Fitting the ST model for black carbon in Denver, CO
Sheena Martenies
01/05/2021
1 Introduction
In this script, we are building the spatiotemporal prediction model for black carbon in the Denver metro area. This model uses data from 5 sampling campaigns (as detailed in 15_Exploring_BC_Data.R). Two models are developed here: one which uses one temporal basis function (Model A) and one that uses two temporal basis functions (Model B).
Updates to this model are in response to reviewer comments on an earlier manuscript version submitted to ES&T.
A summary table at the end of the document summarizes the performance for each of these models. Overall, Model B performs well and allows us to predict BC concentrations at a weekly temporal resolution from 2009 through 2020.
## Loading required package: Matrix2 Developing data sets
2.1 Denver data set
Here I’m reading in the filter data sets and getting them into the right shape for the SpatioTemporal package. Here the ST covariates cover the entire study period (2009 - 2020).
data_name <- "Combined_Filter_Data_AEA.csv"
all_data <- read_csv(here::here("Data", data_name))## Parsed with column specification:
## cols(
## .default = col_double(),
## filter_id = col_character(),
## campaign = col_character(),
## StartDateTimeLocal = col_date(format = ""),
## EndDateTimeLocal = col_date(format = ""),
## sample_week = col_date(format = ""),
## As_ug_m3 = col_logical(),
## Cd_ug_m3 = col_logical(),
## In_ug_m3 = col_logical(),
## Ni_ug_m3 = col_logical(),
## Sn_ug_m3 = col_logical(),
## Te_ug_m3 = col_logical(),
## st_week = col_date(format = ""),
## units_pm = col_character(),
## units_bc = col_character(),
## units_no2 = col_character(),
## units_temp = col_character(),
## units_smoke = col_character(),
## nn3_bc = col_logical(),
## area_bc = col_logical(),
## idw_bc = col_logical()
## )## See spec(...) for full column specifications.#' Select a "calibrated" version of the data
#' For now, go with Deming regression-- accounts for variability in the
#' monitor and the UPAS data and the temporal mismatch in TWAs
all_data$bc_ug_m3 <- all_data$bc_ug_m3_dem
all_data$pm_ug_m3 <- all_data$pm_ug_m3_dem
#' List of sites that failed preliminary screening (see 15_Exploring_BC_Data.R)
#' These are all the sites (by campaign) where BC measurements had a CV > 30%
#' Note that all of these are in Campaign 4
cv_drop <- read_csv(here::here("Data/Dropped_Sites_by_Campaign.csv"))## Parsed with column specification:
## cols(
## site_id = col_double(),
## campaign = col_character(),
## n_filters = col_double(),
## mean = col_double(),
## sd = col_double(),
## cv = col_double(),
## median = col_double(),
## IQR = col_double(),
## range = col_double(),
## cv_gt_50 = col_double(),
## cv_gt_30 = col_double()
## )filter_data <- filter(all_data, is.na(filter_id) | filter_id != "080310027") %>%
filter(is.na(indoor) | indoor == 0) %>%
filter(is.na(is_blank) | is_blank == 0) %>%
#' QA filters
filter(is.na(bc_below_lod) | bc_below_lod == 0) %>%
filter(is.na(negative_pm_mass) | negative_pm_mass == 0) %>%
filter(is.na(potential_contamination) | potential_contamination == 0)
dist_data <- bind_rows(filter(filter_data, is.na(campaign) | campaign != "Campaign4"),
filter(filter_data, campaign == "Campaign4" & !(site_id %in% cv_drop$site_id))) %>%
arrange(st_week)
head(dist_data$sample_week)## [1] NA NA NA NA NA NAtail(dist_data$sample_week)## [1] NA NA NA NA NA NAhead(dist_data$st_week)## [1] "2008-12-29" "2008-12-29" "2008-12-29" "2008-12-29" "2008-12-29"
## [6] "2008-12-29"tail(dist_data$st_week)## [1] "2020-12-28" "2020-12-28" "2020-12-28" "2020-12-28" "2020-12-28"
## [6] "2020-12-28"length(unique(dist_data$filter_id))## [1] 848summary(dist_data$bc_ug_m3)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.51 0.99 1.17 1.27 1.46 2.27 41600quantile(dist_data$bc_ug_m3, probs = c(0.05, 0.95), na.rm = T)## 5% 95%
## 0.8689785 2.0648167sd(dist_data$bc_ug_m3, na.rm = T)## [1] 0.3498975sd(dist_data$bc_ug_m3, na.rm = T) / mean(dist_data$bc_ug_m3, na.rm = T)## [1] 0.2764677#' All of the dropped data should be in Campaign 4
dropped_data <- filter(filter_data, campaign == "Campaign4" & site_id %in% cv_drop$site_id)
length(unique(dropped_data$filter_id))## [1] 97summary(dropped_data$bc_ug_m3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.7752 0.6319 1.4887 1.1935 1.9480 2.1754quantile(dropped_data$bc_ug_m3, probs = c(0.05, 0.95))## 5% 95%
## -0.3137664 2.1225874sd(dropped_data$bc_ug_m3)## [1] 0.9202973sd(dropped_data$bc_ug_m3) / mean(dropped_data$bc_ug_m3)## [1] 0.7711097#' central site data
central_data <- filter(all_data, filter_id == "080310027") %>%
arrange(st_week)
head(central_data$sample_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(central_data$sample_week)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"head(central_data$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(central_data$st_week)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"#' Add prefix to site_ids
unique(dist_data$site_id)## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68dist_data$site_id <- paste0("d_", dist_data$site_id)
unique(dist_data$site_id)## [1] "d_1" "d_2" "d_3" "d_4" "d_5" "d_6" "d_7" "d_8" "d_9" "d_10"
## [11] "d_11" "d_12" "d_13" "d_14" "d_15" "d_16" "d_17" "d_18" "d_19" "d_20"
## [21] "d_21" "d_22" "d_23" "d_24" "d_25" "d_26" "d_27" "d_28" "d_29" "d_30"
## [31] "d_31" "d_32" "d_33" "d_34" "d_35" "d_36" "d_37" "d_38" "d_39" "d_40"
## [41] "d_41" "d_42" "d_43" "d_44" "d_45" "d_46" "d_47" "d_48" "d_49" "d_50"
## [51] "d_51" "d_52" "d_53" "d_54" "d_55" "d_56" "d_57" "d_58" "d_59" "d_60"
## [61] "d_61" "d_62" "d_63" "d_64" "d_65" "d_66" "d_67" "d_68"unique(central_data$site_id)## [1] 16central_data$site_id <- "central"
unique(central_data$site_id)## [1] "central"#' put these data sets back together
lur_data <- bind_rows(dist_data, central_data) %>%
arrange(site_id, st_week)
head(lur_data$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(lur_data$st_week)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"#' Log-transformed BC observations with NA values for the missing dates
#' Note that campaign 5 has duplicate measures, so we need to average them
bc_obs <- lur_data %>%
dplyr::select(site_id, st_week, bc_ug_m3) %>%
group_by(site_id, st_week) %>%
summarize(bc_ug_m3 = mean(bc_ug_m3, na.rm = T)) %>%
mutate(bc_ug_m3 = ifelse(is.nan(bc_ug_m3), NA, bc_ug_m3)) %>%
mutate(log_bc = log(bc_ug_m3)) %>%
dplyr::select(-bc_ug_m3)
#' Pivot to a wide data frame
bc_obs2 <- bc_obs %>%
pivot_wider(id_cols = st_week, names_from = site_id, values_from = log_bc) %>%
# names_from = site_id, values_from = bc_ug_m3) %>%
as.data.frame() %>%
arrange(st_week)
rownames(bc_obs2) <- bc_obs2$st_week
bc_obs2$st_week <- NULL
bc_obs2 <- as.matrix(bc_obs2)
bc_weeks <- rownames(bc_obs2)
tail(bc_weeks)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"#' Spatial covariates (scaled)
load(here::here("Results", "BC_LASSO_Model_5C.Rdata"))
lasso_coef_df <- data.frame(name = log_bc_lasso_coef4@Dimnames[[1]][log_bc_lasso_coef4@i + 1],
coefficient = log_bc_lasso_coef4@x)
covars <- as.character(lasso_coef_df$name[-1])
covars ## [1] "tree_cover_100" "impervious_2500" "low_int_100" "med_int_50"
## [5] "med_int_2500" "high_int_100" "pop_den_50" "dist_m_cafos"
## [9] "dist_m_og_wells" "dist_m_wwtp" "aadt_100" "aadt_250"
## [13] "aadt_1000"covar_fun <- paste("~", paste(covars, collapse = " + "))
covar_fun## [1] "~ tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 + med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos + dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000"bc_sp_cov <- dplyr::select(lur_data, site_id, lon, lat, all_of(covars)) %>%
distinct() %>%
mutate_at(.vars = vars(all_of(covars)),
scale)
bc_sp_cov <- bc_sp_cov %>%
rename(ID = site_id) %>%
mutate(lon2 = lon, lat2 = lat) %>%
st_as_sf(coords = c("lon2", "lat2"), crs = ll_wgs84) %>%
st_transform(crs = albers)
sp_coords <- do.call(rbind, st_geometry(bc_sp_cov)) %>% as_tibble()
names(sp_coords) <- c("x", "y")
bc_sp_cov <- bind_cols(bc_sp_cov, sp_coords) %>%
st_set_geometry(NULL) %>%
as.data.frame()
bc_sp_cov$type <- ifelse(bc_sp_cov$ID == "central", "central", "dist")
bc_sp_cov$type <- as.factor(bc_sp_cov$type)
cor(bc_sp_cov[,covars])## tree_cover_100 impervious_2500 low_int_100 med_int_50
## tree_cover_100 1.000000000 0.05592620 0.227568688 -0.41656599
## impervious_2500 0.055926204 1.00000000 0.046627086 -0.13878441
## low_int_100 0.227568688 0.04662709 1.000000000 -0.11614316
## med_int_50 -0.416565993 -0.13878441 -0.116143165 1.00000000
## med_int_2500 0.005240397 0.89874242 0.172586606 -0.13029847
## high_int_100 -0.335160074 0.48006634 -0.444088483 0.11084350
## pop_den_50 0.201795804 0.39972348 0.483551971 0.05903168
## dist_m_cafos 0.262782544 -0.02217726 -0.161007912 -0.09924381
## dist_m_og_wells 0.304301808 -0.11331048 0.008034836 -0.21284246
## dist_m_wwtp 0.140092353 -0.49234185 -0.075238446 0.02725541
## aadt_100 -0.176079221 0.40911376 -0.255918459 0.06413761
## aadt_250 -0.007357196 0.42742607 -0.261080670 0.03197137
## aadt_1000 0.085749987 0.43050282 -0.112576947 -0.19807436
## med_int_2500 high_int_100 pop_den_50 dist_m_cafos
## tree_cover_100 0.005240397 -0.33516007 0.20179580 0.26278254
## impervious_2500 0.898742424 0.48006634 0.39972348 -0.02217726
## low_int_100 0.172586606 -0.44408848 0.48355197 -0.16100791
## med_int_50 -0.130298466 0.11084350 0.05903168 -0.09924381
## med_int_2500 1.000000000 0.39325577 0.51099015 -0.24164427
## high_int_100 0.393255768 1.00000000 -0.10833432 -0.01905821
## pop_den_50 0.510990152 -0.10833432 1.00000000 -0.11900734
## dist_m_cafos -0.241644265 -0.01905821 -0.11900734 1.00000000
## dist_m_og_wells -0.251552596 -0.21708310 -0.02501575 0.61986048
## dist_m_wwtp -0.456542801 -0.30637142 -0.14260842 0.15352999
## aadt_100 0.261463299 0.71529092 -0.08808862 0.10155682
## aadt_250 0.276833428 0.60158329 -0.10795874 0.23112012
## aadt_1000 0.349181488 0.40408074 0.01155408 0.12765618
## dist_m_og_wells dist_m_wwtp aadt_100 aadt_250
## tree_cover_100 0.304301808 0.14009235 -0.17607922 -0.007357196
## impervious_2500 -0.113310476 -0.49234185 0.40911376 0.427426070
## low_int_100 0.008034836 -0.07523845 -0.25591846 -0.261080670
## med_int_50 -0.212842464 0.02725541 0.06413761 0.031971369
## med_int_2500 -0.251552596 -0.45654280 0.26146330 0.276833428
## high_int_100 -0.217083104 -0.30637142 0.71529092 0.601583294
## pop_den_50 -0.025015754 -0.14260842 -0.08808862 -0.107958743
## dist_m_cafos 0.619860480 0.15352999 0.10155682 0.231120116
## dist_m_og_wells 1.000000000 0.35073325 -0.06090070 0.109442838
## dist_m_wwtp 0.350733252 1.00000000 -0.20875176 -0.024182445
## aadt_100 -0.060900696 -0.20875176 1.00000000 0.806259695
## aadt_250 0.109442838 -0.02418244 0.80625970 1.000000000
## aadt_1000 0.091175296 -0.05838643 0.50962569 0.635217540
## aadt_1000
## tree_cover_100 0.08574999
## impervious_2500 0.43050282
## low_int_100 -0.11257695
## med_int_50 -0.19807436
## med_int_2500 0.34918149
## high_int_100 0.40408074
## pop_den_50 0.01155408
## dist_m_cafos 0.12765618
## dist_m_og_wells 0.09117530
## dist_m_wwtp -0.05838643
## aadt_100 0.50962569
## aadt_250 0.63521754
## aadt_1000 1.00000000#' NO2, and smoke spatiotemporal predictors
#' NO2 at each site estimated using IDW
bc_st_no2 <- dplyr::select(lur_data, site_id, st_week, idw_no2) %>%
distinct() %>%
arrange(st_week) %>%
pivot_wider(id_cols = st_week,
names_from = site_id, values_from = idw_no2) %>%
as.data.frame()
rownames(bc_st_no2) <- bc_st_no2$st_week
bc_st_no2$st_week <- NULL
bc_st_no2 <- as.matrix(bc_st_no2)
no2_weeks <- rownames(bc_st_no2)
tail(no2_weeks)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"setdiff(bc_weeks, no2_weeks)## character(0)#' Smoke day indicator variable based on WFS paper method (see Brey and Fischer, 2016)
#' based on any smoke variable in the area using the 2 sd cut-off (area_smoke_2sd)
bc_st_smk <- dplyr::select(lur_data, site_id, st_week, area_smoke_2sd) %>%
distinct() %>%
arrange(st_week) %>%
pivot_wider(id_cols = st_week,
names_from = site_id, values_from = area_smoke_2sd) %>%
as.data.frame()
rownames(bc_st_smk) <- bc_st_smk$st_week
bc_st_smk$st_week <- NULL
bc_st_smk <- as.matrix(bc_st_smk)
smk_weeks <- rownames(bc_st_smk)
tail(smk_weeks)## [1] "2020-11-23" "2020-11-30" "2020-12-07" "2020-12-14" "2020-12-21"
## [6] "2020-12-28"setdiff(bc_weeks, smk_weeks)## character(0)Create a new denver.data object.
denver.data <- createSTdata(obs = bc_obs2,
covars = bc_sp_cov,
SpatioTemporal = list(bc_st_no2 = bc_st_no2,
bc_st_smk = bc_st_smk))
D <- createDataMatrix(denver.data)
denver.data## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 231 (observed: 231)
## No. obs: 1021
##
## Only constant temporal trend,with dates:
## 2016-02-08 to 2020-11-30
##
## 19 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "med_int_2500"
## [9] "high_int_100" "pop_den_50" "dist_m_cafos" "dist_m_og_wells"
## [13] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"
## [17] "x" "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-062.2 Creating the NO2, BC, and PM data objects
Next we need to combine monitoring data for NO2, BC, and PM2.5 in order to calculate the time trends.
First up: NO2
Plot the time trends for NO2 at the monitors
Original units (ppb)
no2_ts <- ggplot(data = no2_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
no2_ts## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Log-transformed NO2
no2_ts2 <- ggplot(data = no2_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
no2_ts2## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data. Use the log-transformed NO2 here
#' NO2 observations
no2_obs <- no2_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(no2 = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = no2) %>%
as.data.frame() %>%
arrange(st_week)
head(no2_obs)head(no2_obs$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(no2_obs$st_week)## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"# NO2 SP object
no2_data2 <- read_csv(here::here("Data", "Monitor_NO2_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_no2")) %>%
filter(monitor_id %in% colnames(no2_obs))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.no2_coords <- do.call(rbind, st_geometry(no2_data2)) %>% as_tibble()
names(no2_coords) <- c("x", "y")
no2_sp_lonlat <- no2_data2 %>%
st_transform(crs = ll_wgs84)
no2_coords2 <- do.call(rbind, st_geometry(no2_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
no2_sp_cov <- st_set_geometry(no2_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
no2_sp_cov <- bind_cols(no2_sp_cov, no2_coords, no2_coords2) %>%
as.data.frame() %>%
distinct()Next, BC at the central site monitor
#' BC central site concentrations
bc_cent_data <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
filter(!is.na(Arithmetic_Mean)) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(County_Code = str_sub(monitor_id, start = 3, end = 5)) %>%
filter(County_Code %in% counties) %>%
arrange(Date_Local, monitor_id) %>%
mutate(Arithmetic_Mean = ifelse(Arithmetic_Mean <= 0.005, NA, Arithmetic_Mean))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Pollutant_Standard = col_logical(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_logical(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.head(bc_cent_data$Date_Local)## [1] "2016-02-08" "2016-02-09" "2016-02-10" "2016-02-11" "2016-02-12"
## [6] "2016-02-13"tail(bc_cent_data$Date_Local)## [1] "2020-11-25" "2020-11-26" "2020-11-27" "2020-11-28" "2020-11-29"
## [6] "2020-11-30"#' Add an identifier to differentiate these measurements from collocated pollutants at the same site
bc_cent_data$monitor_id <- paste0(bc_cent_data$monitor_id, "_bc")
unique(bc_cent_data$monitor_id)## [1] "080310027_bc"Plot the time trends for BC at the central site monitor
Original units (ug/m3)
bc_ts <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
#facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
bc_ts## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Log-transformed BC
bc_ts2 <- ggplot(data = bc_cent_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
#facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
bc_ts2## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data
Use log-transformed data here
#' BC observations
bc_cent_obs <- bc_cent_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(bc = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = bc) %>%
as.data.frame() %>%
arrange(st_week)
head(bc_cent_obs)head(bc_cent_obs$st_week)## [1] "2016-02-08" "2016-02-15" "2016-02-22" "2016-02-29" "2016-03-07"
## [6] "2016-03-14"tail(bc_cent_obs$st_week)## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"# BC Central Site SP object
bc_cent_data2 <- read_csv(here::here("Data", "Monitor_BC_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_bc")) %>%
filter(monitor_id %in% colnames(bc_cent_obs))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Pollutant_Standard = col_logical(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_logical(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.bc_cent_coords <- do.call(rbind, st_geometry(bc_cent_data2)) %>% as_tibble()
names(bc_cent_coords) <- c("x", "y")
bc_cent_sp_lonlat <- bc_cent_data2 %>%
st_transform(crs = ll_wgs84)
bc_cent_coords2 <- do.call(rbind, st_geometry(bc_cent_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
bc_cent_sp_cov <- st_set_geometry(bc_cent_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
bc_cent_sp_cov <- bind_cols(bc_cent_sp_cov, bc_cent_coords, bc_cent_coords2) %>%
as.data.frame() %>%
distinct()Lastly, PM2.5 from select monitors with long-term records
Plot the time trends for PM2.5 at the monitors
Original units (ug/m3)
pm_ts <- ggplot(data = pm_data, aes(x = Date_Local, y = Arithmetic_Mean)) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
pm_ts## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Log-transformed PM2.5
pm_ts2 <- ggplot(data = pm_data, aes(x = Date_Local, y = log(Arithmetic_Mean))) +
geom_point(aes(color = as.factor(monitor_id)), shape = 20, size = 1) +
geom_smooth(se = F, color = "black", method = "gam") +
scale_color_viridis(name = "Monitor ID", discrete = T) +
scale_x_date(date_breaks = "6 months", date_labels = "%m-%Y") +
facet_wrap(. ~ monitor_id, scales = "fixed", ncol = 2) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
simple_theme
pm_ts2## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
Organize the observations and spatial data
Use the log-transformed data here
#' PM2.5 observations
pm_obs <- pm_data %>%
st_set_geometry(NULL) %>%
dplyr::select(monitor_id, Date_Local, Arithmetic_Mean) %>%
mutate(st_week = as.character(as.Date(cut(as.Date(Date_Local), "week")))) %>%
filter(st_week %in% bc_weeks) %>%
group_by(monitor_id, st_week) %>%
mutate(log_mean = log(Arithmetic_Mean)) %>%
summarize(pm = mean(log_mean, na.rm=T)) %>%
pivot_wider(id_cols = st_week,
names_from = monitor_id, values_from = pm) %>%
as.data.frame() %>%
arrange(st_week)
head(pm_obs)head(pm_obs$st_week)## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(pm_obs$st_week)## [1] "2020-10-26" "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23"
## [6] "2020-11-30"# PM SP object
pm_data2 <- read_csv(here::here("Data", "Monitor_PM_Data_AEA.csv")) %>%
mutate(year = year(Date_Local)) %>%
dplyr::select(-year) %>%
st_as_sf(wkt = "WKT", crs = albers) %>%
mutate(monitor_id = paste0(monitor_id, "_pm")) %>%
filter(monitor_id %in% colnames(pm_obs))## Parsed with column specification:
## cols(
## .default = col_character(),
## Parameter_Code = col_double(),
## POC = col_double(),
## Date_Local = col_date(format = ""),
## Observation_Count = col_double(),
## Observation_Percent = col_double(),
## Arithmetic_Mean = col_double(),
## Max_Value = col_double(),
## `1st_Max_Hour` = col_double(),
## AQI = col_double(),
## Method_Code = col_double(),
## Date_of_Last_Change = col_date(format = "")
## )## See spec(...) for full column specifications.pm_coords <- do.call(rbind, st_geometry(pm_data2)) %>% as_tibble()
names(pm_coords) <- c("x", "y")
pm_sp_lonlat <- pm_data2 %>%
st_transform(crs = ll_wgs84)
pm_coords2 <- do.call(rbind, st_geometry(pm_sp_lonlat)) %>%
as_tibble() %>% setNames(c("lon","lat"))
pm_sp_cov <- st_set_geometry(pm_data2, NULL) %>%
dplyr::select(monitor_id) %>%
rename(ID = monitor_id)
pm_sp_cov <- bind_cols(pm_sp_cov, pm_coords, pm_coords2) %>%
as.data.frame() %>%
distinct()Combine the observations and the spatial covariates Scale the pollutant measurements
#' Join log-transformed observations
pol_obs <- left_join(no2_obs, bc_cent_obs, by = "st_week") %>%
left_join(pm_obs, by = "st_week")
glimpse(pol_obs)## Rows: 622
## Columns: 18
## $ st_week <chr> "2008-12-29", "2009-01-05", "2009-01-12", "2009-01-19…
## $ `080013001_no2` <dbl> 2.4313422, 2.3337153, 3.0953820, 2.4960668, 2.6445846…
## $ `080310002_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310026_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_no2` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_bc` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080010006_pm` <dbl> 1.394118, 1.247789, 1.934844, 2.298456, 1.598321, 1.7…
## $ `080010008_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080050005_pm` <dbl> 1.2490759, 0.6817687, 1.3697744, 1.7730816, 1.1939225…
## $ `080130003_pm` <dbl> 0.9679299, 0.8243293, 2.1047287, 1.7212142, 1.9059908…
## $ `080130012_pm` <dbl> 1.0919008, 0.8908546, 1.2511276, 1.2703657, 0.8047190…
## $ `080310002_pm` <dbl> 1.335209, 1.188381, 1.767997, 2.027844, 1.486213, 1.7…
## $ `080310023_pm` <dbl> 1.661197, 1.513911, 1.919262, 1.996528, 1.909773, 1.7…
## $ `080310025_pm` <dbl> 1.1537863, 0.9359011, 1.3262687, 1.9969858, 1.6521301…
## $ `080310026_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310027_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ `080310028_pm` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…rownames(pol_obs) <- pol_obs$st_week
pol_obs$st_week <- NULL
pol_obs <- as.matrix(pol_obs)
#' Scale the measurements to avoid issues where units are different
#' Remember! These have been log-transformed before scaling
#' Scale the data by columns (monitors)
rnames <- rownames(pol_obs)
cnames <- colnames(pol_obs)
pol_obs <- apply(pol_obs, 2, scale)
rownames(pol_obs) <- rnames
colnames(pol_obs) <- cnames
head(rownames(pol_obs))## [1] "2008-12-29" "2009-01-05" "2009-01-12" "2009-01-19" "2009-01-26"
## [6] "2009-02-02"tail(rownames(pol_obs))## [1] "2020-11-02" "2020-11-09" "2020-11-16" "2020-11-23" "2020-11-30"
## [6] "2020-12-07"head(colnames(pol_obs))## [1] "080013001_no2" "080310002_no2" "080310026_no2" "080310027_no2"
## [5] "080310028_no2" "080310027_bc"class(pol_obs)## [1] "matrix" "array"dim(pol_obs)## [1] 622 17#' SP variables
pol_sp_cov <- bind_rows(no2_sp_cov, bc_cent_sp_cov, pm_sp_cov)
pol_sp_covCreate the data object
pol_STdata <- createSTdata(obs = pol_obs,
covars = pol_sp_cov)
print(pol_STdata)## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Only constant temporal trend,with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.2.3 Generate the temporal trend
The cross validation results for generating the time trends suggest 2 basis functions is the way to go, but let’s see if that holds up with our data
D_pol <- createDataMatrix(pol_STdata)
#D_pol
n_years <- length(2009:2020)
n_years*4## [1] 48n_years*8## [1] 96#' Trying 4 df per year
pol.SVD.cv.4py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 4)
print(pol.SVD.cv.4py)## Result of SVDsmoothCV, average of CV-statistics:
## MSE R2 AIC BIC
## n.basis.0 0.9975610 0.0000000 0.9987781 5.095695
## n.basis.1 0.8013682 0.1967023 -96.6124965 -88.418663
## n.basis.2 0.6358644 0.3625750 -225.3897791 -213.099029
## n.basis.3 0.6297167 0.3687387 -229.4297802 -213.042114
## n.basis.4 0.6084604 0.3900698 -240.5027300 -220.018147plot(pol.SVD.cv.4py)
#' Trying 8 df per year
pol.SVD.cv.8py <- SVDsmoothCV(D_pol, 0:4, df = n_years * 8)
print(pol.SVD.cv.8py)## Result of SVDsmoothCV, average of CV-statistics:
## MSE R2 AIC BIC
## n.basis.0 0.9975610 0.0000000 0.9987781 5.095695
## n.basis.1 0.7630671 0.2350911 -120.2211615 -112.027328
## n.basis.2 0.5825979 0.4159654 -265.3105795 -253.019830
## n.basis.3 0.5739736 0.4246129 -271.3027690 -254.915103
## n.basis.4 0.5615078 0.4371173 -278.8000243 -258.315442plot(pol.SVD.cv.8py)
2.4 Plot the time trends
I tried using 1 and 2 basis functions for the time trends. I also checked out 4 and 8 df per year
2.4.1 1 basis function
First, see what the time trends for the monitoring data look like with 1 basis function.
4 df per year
pol_STdata1 <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 4)## Replacing existing trend.pol_STdata1## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata1$trend)tail(pol_STdata1$trend)plot(pol_STdata1$trend$date, pol_STdata1$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1$trend$date, pol_STdata1$atrend$V1, col = 1)
8 df per year
pol_STdata1a <- updateTrend(pol_STdata, n.basis = 1, df = n_years * 8)## Replacing existing trend.pol_STdata1a## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata1a$trend)tail(pol_STdata1a$trend)plot(pol_STdata1a$trend$date, pol_STdata1a$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
lines(pol_STdata1a$trend$date, pol_STdata1a$atrend$V1, col = 1)
2.4.2 2 basis functions
Next, let’s see what 2 basis functions look like:
4 df per year
pol_STdata2 <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 4)## Replacing existing trend.pol_STdata2## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata2$trend)tail(pol_STdata2$trend)plot(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V1, col = 1)
lines(pol_STdata2$trend$date, pol_STdata2$trend$V2, col = 2)
8 df per year
pol_STdata2a <- updateTrend(pol_STdata, n.basis = 2, df = n_years * 8)## Replacing existing trend.pol_STdata2a## STdata-object with:
## No. locations: 17 (observed: 16)
## No. time points: 622 (observed: 622)
## No. obs: 6560
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 5 covariate(s):
## [1] "ID" "x" "y" "lon" "lat"
##
## No spatio-temporal covariates.head(pol_STdata2a$trend)tail(pol_STdata2a$trend)plot(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1, pch = 16, cex = 0.5,
xlab = "Date", ylab = "Pollutants (scaled)")
points(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2, pch = 16, cex = 0.5)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V1, col = 1)
lines(pol_STdata2a$trend$date, pol_STdata2a$trend$V2, col = 2)
2.5 Update the Denver data set
The next step is to update the trend data for the denver.data object and see if these trend data fit our observations.
Comparing the plots using one basis function vs. two basis functions, it looks like using two basis functions might be better than 1. There was some residual autocorrelation in the central site monitor data when using just one basis function (see the fourth plot in the first series of plots below). Using 8 degrees of freedom did not help.
Using two basis functions helped somewhat to reduce this autocorrelation, but it still persists. There are still some noticable trends in the residuals for the central monitoring site
2.5.1 1 basis function, 4 df per year
First update the trend data, then plot the trends for three distributed sites and the central BC monitoring site.
denver.data2.1 <- denver.data
denver.data2.1$trend <- pol_STdata1$trend
print(denver.data2.1)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 19 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "med_int_2500"
## [9] "high_int_100" "pop_den_50" "dist_m_cafos" "dist_m_og_wells"
## [13] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"
## [17] "x" "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.1, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_1")
plot(denver.data2.1, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_10")
plot(denver.data2.1, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_43")
plot(denver.data2.1, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.1, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="d_53")
plot(denver.data2.1, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.1, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1, "acf", ID="central")
plot(denver.data2.1, "pacf", ID="central")
2.5.2 1 basis function, 8 df per year
denver.data2.1a <- denver.data
denver.data2.1a$trend <- pol_STdata1a$trend
print(denver.data2.1a)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 19 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "med_int_2500"
## [9] "high_int_100" "pop_den_50" "dist_m_cafos" "dist_m_og_wells"
## [13] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"
## [17] "x" "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.1a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_1")
plot(denver.data2.1a, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_10")
plot(denver.data2.1a, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.1a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_43")
plot(denver.data2.1a, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.1a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="d_53")
plot(denver.data2.1a, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.1a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.1a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.1a, "acf", ID="central")
plot(denver.data2.1a, "pacf", ID="central")
2.5.3 2 basis functions, 4 df per year
denver.data2.2 <- denver.data
denver.data2.2$trend <- pol_STdata2$trend
print(denver.data2.2)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 19 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "med_int_2500"
## [9] "high_int_100" "pop_den_50" "dist_m_cafos" "dist_m_og_wells"
## [13] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"
## [17] "x" "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.2, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_1")
plot(denver.data2.2, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_10")
plot(denver.data2.2, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_10")
plot(denver.data2.2, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.2, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_43")
plot(denver.data2.2, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.2, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="d_53")
plot(denver.data2.2, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.2, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2, "acf", ID="central")
plot(denver.data2.2, "pacf", ID="central")
2.5.4 2 basis functions, 8 df per year
denver.data2.2a <- denver.data
denver.data2.2a$trend <- pol_STdata2a$trend
print(denver.data2.2a)## STdata-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## 19 covariate(s):
## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "med_int_2500"
## [9] "high_int_100" "pop_den_50" "dist_m_cafos" "dist_m_og_wells"
## [13] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"
## [17] "x" "y" "type"
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-06par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_1", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_1")
plot(denver.data2.2a, "res", ID="d_1", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_1")
plot(denver.data2.2a, "pacf", ID="d_1")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_10", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_10")
plot(denver.data2.2a, "res", ID="d_10", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_10")
plot(denver.data2.2a, "pacf", ID="d_10")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_43", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_43")
plot(denver.data2.2a, "res", ID="d_43", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_43")
plot(denver.data2.2a, "pacf", ID="d_43")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="d_53", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend d_53")
plot(denver.data2.2a, "res", ID="d_53", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="d_53")
plot(denver.data2.2a, "pacf", ID="d_53")
par(mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow=TRUE))
plot(denver.data2.2a, "obs", ID="central", xlab="", ylab="BC (log ug/m3)",
main="Temporal trend central")
plot(denver.data2.2a, "res", ID="central", xlab="", ylab="BC (log ug/m3)")
plot(denver.data2.2a, "acf", ID="central")
plot(denver.data2.2a, "pacf", ID="central")
3 Testing possible models
Building on the results of the previous modeling efforts, I’m going to update “Model 4.2”, which used an ’exp covariance structure for the error term and iid for the beta fields.
3.1 Model A: 1 temporal trend (basis) function (4 df per year)
3.1.1 Create the model object
For this version of the model, use iid for cov.beta (beta, beta1) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.
denver.data.A <- denver.data2.1
names(denver.data.A$covars)## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "med_int_2500"
## [9] "high_int_100" "pop_den_50" "dist_m_cafos" "dist_m_og_wells"
## [13] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"
## [17] "x" "y" "type"LUR.A <- list(covar_fun, covar_fun)
cov.beta.A <- list(covf = c("iid", "iid"), nugget = T)
cov.nu.A <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.A <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")
denver.model.A <- createSTmodel(denver.data.A, LUR = LUR.A,
ST = c("bc_st_no2", "bc_st_smk"),
cov.beta = cov.beta.A, cov.nu = cov.nu.A,
locations = locations.A)## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.denver.model.A## STmodel-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 1 basis function(s):
## [1] "V1"
## with dates:
## 2008-12-29 to 2020-12-07
##
## Models for the beta-fields are:
## $const
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos +
## dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## $V1
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos +
## dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## Covariance model for the beta-field(s):
## Covariance type(s): iid, iid
## Nugget: Yes, Yes
## Covariance model for the nu-field(s):
## Covariance type: exp
## Nugget: ~1
## Random effect: No
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-063.1.2 Estimate model parameters
- nugget values can range from -1 to -6
- range values can range from 0 to 4
- try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
- make sure the starting values are pretty different
set.seed(1000)
names <- loglikeSTnames(denver.model.A, all=FALSE)
names## [1] "log.nugget.const.iid" "log.nugget.V1.iid"
## [3] "nu.log.range.exp" "nu.log.sill.exp"
## [5] "nu.log.nugget.(Intercept).exp"# x.init.A <- cbind(c(0, 0, 0, 0, 0),
# c(-1, -1, 0, -1, -1),
# c(-1, -1, 0, -5, -1),
# c(-5, -5, 0, -1, -5),
# c(-5, -5, 0, -5, -5),
# c(-1, -1, 2, -1, -1),
# c(-1, -1, 2, -5, -1),
# c(-5, -5, 2, -1, -5),
# c(-5, -5, 2, -5, -5),
# c(-1, -1, 4, -1, -1),
# c(-1, -1, 4, -5, -1),
# c(-5, -5, 4, -1, -5),
# c(-5, -5, 4, -5, -5))
x.init.A <- cbind(c(0, 0, 0, 0, 0),
c(-1, -1, 4, -5, -1),
c(-5, -5, 4, -1, -5),
c(-1, -1, 8, -5, -1),
c(-5, -5, 8, -1, -5),
c(-1, -1, 12, -5, -1),
c(-5, -5, 12, -1, -5),
c(-7, -7, 12, -3, -3))
rownames(x.init.A) <- loglikeSTnames(denver.model.A, all=FALSE)
x.init.A## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## log.nugget.const.iid 0 -1 -5 -1 -5 -1 -5 -7
## log.nugget.V1.iid 0 -1 -5 -1 -5 -1 -5 -7
## nu.log.range.exp 0 4 4 8 8 12 12 12
## nu.log.sill.exp 0 -5 -1 -5 -1 -5 -1 -3
## nu.log.nugget.(Intercept).exp 0 -1 -5 -1 -5 -1 -5 -3#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.A <- estimate.STmodel(denver.model.A, x.init.A)## Optimisation using starting value 1/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= 459.08 |proj g|= 15
## At iterate 10 f = -1309.8 |proj g|= 1.2358
## At iterate 20 f = -1310 |proj g|= 0.0064564
##
## iterations 21
## function evaluations 35
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 2
## norm of the final projected gradient 0.00646143
## final function value -1310.01
##
## F = -1310.01
## final value -1310.007198
## converged## Optimisation using starting value 2/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -331.87 |proj g|= 14
## At iterate 10 f = -1306.2 |proj g|= 1.925
## ys=-2.142e-01 -gs= 4.291e-02, BFGS update SKIPPED
## At iterate 20 f = -1307.4 |proj g|= 11.314
## At iterate 30 f = -1310 |proj g|= 0.5797
## At iterate 40 f = -1402.8 |proj g|= 20.928
## At iterate 50 f = -1572.4 |proj g|= 0.014546
##
## iterations 51
## function evaluations 93
## segments explored during Cauchy searches 54
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.014546
## final function value -1572.4
##
## F = -1572.4
## final value -1572.396317
## converged## Optimisation using starting value 3/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -458.66 |proj g|= 14
## At iterate 10 f = -1304.7 |proj g|= 1.9297
## At iterate 20 f = -1308.9 |proj g|= 10.68
## ys=-2.936e+00 -gs= 9.218e-01, BFGS update SKIPPED
## At iterate 30 f = -1540 |proj g|= 20.292
##
## iterations 39
## function evaluations 53
## segments explored during Cauchy searches 40
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00692961
## final function value -1572.4
##
## F = -1572.4
## final value -1572.397522
## converged## Optimisation using starting value 4/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -332.17 |proj g|= 14
## At iterate 10 f = -1571.2 |proj g|= 3.1726
## At iterate 20 f = -1576.3 |proj g|= 0.32053
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 29
## function evaluations 79
## segments explored during Cauchy searches 33
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00351314
## final function value -1576.3
##
## F = -1576.3
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1576.303657
## converged## Optimisation using starting value 5/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -546.73 |proj g|= 14
## At iterate 10 f = -1576 |proj g|= 1.1049
## At iterate 20 f = -1576.3 |proj g|= 0.08239
##
## iterations 24
## function evaluations 46
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0133713
## final function value -1576.3
##
## F = -1576.3
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1576.303693
## converged## Optimisation using starting value 6/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -335.65 |proj g|= 14
## At iterate 10 f = -1576.1 |proj g|= 6.8273
##
## iterations 19
## function evaluations 37
## segments explored during Cauchy searches 22
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.076101
## final function value -1576.3
##
## F = -1576.3
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1576.303485
## converged## Optimisation using starting value 7/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1296.8 |proj g|= 14
## At iterate 10 f = -1576.3 |proj g|= 0.36275
## At iterate 20 f = -1576.3 |proj g|= 0.06807
##
## iterations 26
## function evaluations 35
## segments explored during Cauchy searches 30
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0533678
## final function value -1576.3
##
## F = -1576.3
## final value -1576.303711
## converged## Optimisation using starting value 8/8## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1228.8 |proj g|= 12
## At iterate 10 f = -1576.3 |proj g|= 1.8747
## At iterate 20 f = -1576.3 |proj g|= 0.082774
##
## iterations 22
## function evaluations 28
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0441012
## final function value -1576.3
##
## F = -1576.3
## final value -1576.303691
## convergedprint(est.denver.model.A)## Optimisation for STmodel with 8 starting points.
## Results: 5 converged, 3 not converged, 0 failed.
## Best result for starting point 7, optimisation has converged
##
## No fixed parameters.
##
## Estimated parameters for all starting point(s):
## [,1] [,2] [,3]
## gamma.bc_st_no2 0.018662027683 0.011947334 0.011948431
## gamma.bc_st_smk -0.045770923734 -0.024770380 -0.024768403
## alpha.const.(Intercept) -0.144213749735 0.043061298 0.043040010
## alpha.const.tree_cover_100 -0.001271761777 0.003719008 0.003718359
## alpha.const.impervious_2500 0.052029988535 0.044170288 0.044168348
## alpha.const.low_int_100 -0.001653106787 0.009016166 0.009015446
## alpha.const.med_int_50 -0.002288302995 0.002734351 0.002734279
## alpha.const.med_int_2500 -0.010397390358 -0.017786781 -0.017786385
## alpha.const.high_int_100 0.010312389425 0.016163245 0.016162803
## alpha.const.pop_den_50 -0.015664209451 -0.013249908 -0.013248458
## alpha.const.dist_m_cafos -0.045208992335 -0.037819717 -0.037820572
## alpha.const.dist_m_og_wells 0.010367781069 0.002794518 0.002795706
## alpha.const.dist_m_wwtp 0.005191651255 0.006739532 0.006738714
## alpha.const.aadt_100 0.034921854611 0.034342561 0.034342776
## alpha.const.aadt_250 0.008441674556 0.006778407 0.006777233
## alpha.const.aadt_1000 -0.001209607701 0.001401806 0.001401659
## alpha.V1.(Intercept) 0.186004760519 0.222990414 0.222990813
## alpha.V1.tree_cover_100 -0.020164715529 -0.016421057 -0.016422865
## alpha.V1.impervious_2500 -0.020298784866 -0.038355053 -0.038358586
## alpha.V1.low_int_100 -0.016748233869 -0.008141071 -0.008142264
## alpha.V1.med_int_50 -0.022951471286 -0.017767597 -0.017767805
## alpha.V1.med_int_2500 0.026917676047 0.025099230 0.025101869
## alpha.V1.high_int_100 0.004004735689 0.009705730 0.009702999
## alpha.V1.pop_den_50 -0.000008909438 0.002641688 0.002640935
## alpha.V1.dist_m_cafos 0.011071875338 0.007011703 0.007010877
## alpha.V1.dist_m_og_wells -0.004784189400 0.001697012 0.001696335
## alpha.V1.dist_m_wwtp -0.009890449702 -0.011769544 -0.011770430
## alpha.V1.aadt_100 -0.046053857986 -0.041872920 -0.041868833
## alpha.V1.aadt_250 0.009842151159 0.007113844 0.007113932
## alpha.V1.aadt_1000 -0.007732179733 -0.005245424 -0.005245087
## log.nugget.const.iid -15.000000000000 -14.999371995 -14.998712098
## log.nugget.V1.iid -15.000000000000 -14.999810884 -14.463366983
## nu.log.range.exp 0.000000000000 12.464310836 12.463721125
## nu.log.sill.exp -4.390684341790 -3.022719588 -3.022834153
## nu.log.nugget.(Intercept).exp -4.135497997773 -4.771806577 -4.771972935
## [,4] [,5] [,6]
## gamma.bc_st_no2 0.012195145 0.012191546 0.012201358
## gamma.bc_st_smk -0.023368096 -0.023376519 -0.023355592
## alpha.const.(Intercept) 0.034510427 0.034586281 0.034389759
## alpha.const.tree_cover_100 0.002813683 0.002817342 0.002809626
## alpha.const.impervious_2500 0.044148457 0.044155063 0.044138336
## alpha.const.low_int_100 0.008272737 0.008276440 0.008268586
## alpha.const.med_int_50 0.001462357 0.001462844 0.001465164
## alpha.const.med_int_2500 -0.019261212 -0.019261589 -0.019257039
## alpha.const.high_int_100 0.014684016 0.014687448 0.014683691
## alpha.const.pop_den_50 -0.011372682 -0.011379302 -0.011366384
## alpha.const.dist_m_cafos -0.038056720 -0.038053453 -0.038064354
## alpha.const.dist_m_og_wells 0.002481102 0.002476676 0.002492421
## alpha.const.dist_m_wwtp 0.005861502 0.005865456 0.005856464
## alpha.const.aadt_100 0.033664281 0.033663797 0.033666382
## alpha.const.aadt_250 0.007104613 0.007108579 0.007095515
## alpha.const.aadt_1000 0.002106792 0.002106273 0.002102821
## alpha.V1.(Intercept) 0.226073249 0.226100754 0.226019975
## alpha.V1.tree_cover_100 -0.016992745 -0.016986819 -0.017002648
## alpha.V1.impervious_2500 -0.036778794 -0.036765912 -0.036806336
## alpha.V1.low_int_100 -0.008874293 -0.008869243 -0.008882278
## alpha.V1.med_int_50 -0.016196431 -0.016196237 -0.016203423
## alpha.V1.med_int_2500 0.024306189 0.024301028 0.024314867
## alpha.V1.high_int_100 0.009012393 0.009021579 0.008997412
## alpha.V1.pop_den_50 0.002208200 0.002208468 0.002210902
## alpha.V1.dist_m_cafos 0.006169577 0.006170451 0.006173326
## alpha.V1.dist_m_og_wells 0.001844731 0.001849461 0.001832932
## alpha.V1.dist_m_wwtp -0.010408832 -0.010408312 -0.010415084
## alpha.V1.aadt_100 -0.040115551 -0.040129281 -0.040090907
## alpha.V1.aadt_250 0.007671507 0.007668694 0.007672913
## alpha.V1.aadt_1000 -0.005173673 -0.005174620 -0.005171007
## log.nugget.const.iid -7.797159485 -7.798415070 -7.800270351
## log.nugget.V1.iid -7.521046376 -7.525829277 -7.515414412
## nu.log.range.exp 12.444970355 12.446646679 12.441842084
## nu.log.sill.exp -3.008558413 -3.008570123 -3.008497539
## nu.log.nugget.(Intercept).exp -4.861112033 -4.860671802 -4.861403855
## [,7] [,8]
## gamma.bc_st_no2 0.012187315 0.012191285
## gamma.bc_st_smk -0.023383491 -0.023379810
## alpha.const.(Intercept) 0.034674416 0.034601251
## alpha.const.tree_cover_100 0.002819106 0.002814760
## alpha.const.impervious_2500 0.044164671 0.044151847
## alpha.const.low_int_100 0.008279085 0.008274175
## alpha.const.med_int_50 0.001460604 0.001462762
## alpha.const.med_int_2500 -0.019266096 -0.019261523
## alpha.const.high_int_100 0.014687926 0.014685727
## alpha.const.pop_den_50 -0.011382275 -0.011375796
## alpha.const.dist_m_cafos -0.038050220 -0.038055208
## alpha.const.dist_m_og_wells 0.002471527 0.002480366
## alpha.const.dist_m_wwtp 0.005867949 0.005862805
## alpha.const.aadt_100 0.033662812 0.033665304
## alpha.const.aadt_250 0.007114408 0.007105814
## alpha.const.aadt_1000 0.002106229 0.002104928
## alpha.V1.(Intercept) 0.226138646 0.226089325
## alpha.V1.tree_cover_100 -0.016979945 -0.016990544
## alpha.V1.impervious_2500 -0.036741674 -0.036773466
## alpha.V1.low_int_100 -0.008864162 -0.008872411
## alpha.V1.med_int_50 -0.016194695 -0.016197912
## alpha.V1.med_int_2500 0.024284973 0.024299875
## alpha.V1.high_int_100 0.009035275 0.009017501
## alpha.V1.pop_den_50 0.002211442 0.002210976
## alpha.V1.dist_m_cafos 0.006174037 0.006172359
## alpha.V1.dist_m_og_wells 0.001851383 0.001843365
## alpha.V1.dist_m_wwtp -0.010402080 -0.010408105
## alpha.V1.aadt_100 -0.040142669 -0.040117310
## alpha.V1.aadt_250 0.007667448 0.007669650
## alpha.V1.aadt_1000 -0.005177372 -0.005173688
## log.nugget.const.iid -7.797036706 -7.798340467
## log.nugget.V1.iid -7.526945504 -7.521553058
## nu.log.range.exp 12.449157106 12.446308860
## nu.log.sill.exp -3.008551080 -3.007960669
## nu.log.nugget.(Intercept).exp -4.860112385 -4.861008888
##
## Function value(s):
## [1] 1310.007 1572.396 1572.398 1576.304 1576.304 1576.303 1576.304 1576.3043.1.3 Cross-validation
Define the CV groups
set.seed(1000)
site_idsA <- colnames(denver.model.A$D.beta)[which(str_detect(colnames(denver.model.A$D.beta), "d_") == T)]
unique(colnames(bc_obs2))## [1] "central" "d_1" "d_10" "d_11" "d_12" "d_13" "d_14"
## [8] "d_15" "d_16" "d_17" "d_18" "d_19" "d_2" "d_20"
## [15] "d_21" "d_22" "d_23" "d_24" "d_25" "d_26" "d_27"
## [22] "d_28" "d_29" "d_3" "d_30" "d_31" "d_32" "d_33"
## [29] "d_34" "d_35" "d_36" "d_37" "d_38" "d_39" "d_4"
## [36] "d_40" "d_41" "d_42" "d_43" "d_44" "d_45" "d_46"
## [43] "d_47" "d_48" "d_49" "d_5" "d_50" "d_51" "d_52"
## [50] "d_53" "d_54" "d_55" "d_56" "d_57" "d_58" "d_59"
## [57] "d_6" "d_60" "d_61" "d_62" "d_63" "d_64" "d_65"
## [64] "d_66" "d_67" "d_68" "d_7" "d_8" "d_9"Ind.cv.A <- createCV(denver.model.A, groups = 10, min.dist = .1,
subset = site_idsA)
ID.cv.A <- sapply(split(denver.model.A$obs$ID, Ind.cv.A), unique)
print(sapply(ID.cv.A, length))## 0 1 2 3 4 5 6 7 8 9 10
## 1 7 7 7 6 6 6 6 6 6 6table(Ind.cv.A)## Ind.cv.A
## 0 1 2 3 4 5 6 7 8 9 10
## 231 75 80 97 73 64 93 80 81 66 81I.col.A <- apply(sapply(ID.cv.A, function(x) denver.model.A$locations$ID%in% x), 1,
function(x) if(sum(x)==1) which(x) else 0)
names(I.col.A) <- denver.model.A$locations$ID
print(I.col.A)## central d_1 d_10 d_11 d_12 d_13 d_14 d_15 d_16 d_17
## 1 6 2 3 4 4 4 9 4 10
## d_18 d_19 d_2 d_20 d_21 d_22 d_23 d_25 d_26 d_28
## 8 5 11 3 2 11 4 9 11 9
## d_29 d_3 d_31 d_32 d_33 d_34 d_35 d_36 d_37 d_38
## 2 4 6 8 5 5 10 8 3 8
## d_39 d_4 d_40 d_41 d_42 d_43 d_44 d_45 d_46 d_47
## 6 7 7 3 7 8 11 3 9 5
## d_48 d_49 d_5 d_50 d_51 d_52 d_53 d_54 d_55 d_56
## 9 4 6 9 10 3 7 5 7 6
## d_57 d_58 d_6 d_60 d_61 d_62 d_63 d_64 d_65 d_66
## 2 2 11 5 2 11 10 3 8 6
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 7 10 2 10 0 0 0 0 0par(mfrow=c(1,1))
plot(denver.model.A$locations$long,
denver.model.A$locations$lat,
pch=23+floor(I.col.A/max(I.col.A)+.5), bg=I.col.A,
xlab="Longitude", ylab="Latitude")
map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)## [[1]]
## NULLID starting conditions using previous model without CV:
x.init.A.cv <- coef(est.denver.model.A, pars="cov")[,c("par","init")]
x.init.A.cvRun the model with cross validation
est.denver.A.cv <- estimateCV(denver.model.A, x.init.A.cv, Ind.cv.A)##
## ***************************## Estimation of CV-set 1/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1430.9 |proj g|= 11.381
## At iterate 10 f = -1431.7 |proj g|= 0.27352
##
## iterations 17
## function evaluations 30
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0379531
## final function value -1431.71
##
## F = -1431.71
## final value -1431.711037
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1163.8 |proj g|= 14
## At iterate 10 f = -1431.5 |proj g|= 0.61862
## At iterate 20 f = -1431.7 |proj g|= 0.13907
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 25
## function evaluations 61
## segments explored during Cauchy searches 30
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0694925
## final function value -1431.71
##
## F = -1431.71
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1431.711043
## converged## ##
## ***************************## Estimation of CV-set 2/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1435.9 |proj g|= 1.0766
## At iterate 10 f = -1435.9 |proj g|= 0.020506
##
## iterations 10
## function evaluations 28
## segments explored during Cauchy searches 10
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0205057
## final function value -1435.92
##
## F = -1435.92
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1435.924803
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1169.1 |proj g|= 14
## At iterate 10 f = -1435.9 |proj g|= 0.22268
##
## iterations 19
## function evaluations 34
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00810512
## final function value -1435.92
##
## F = -1435.92
## final value -1435.924810
## converged## ##
## ***************************## Estimation of CV-set 3/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1419.9 |proj g|= 10.14
## At iterate 10 f = -1420.6 |proj g|= 0.23026
##
## iterations 17
## function evaluations 31
## segments explored during Cauchy searches 19
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0225161
## final function value -1420.64
##
## F = -1420.64
## final value -1420.637129
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1149.9 |proj g|= 14
## At iterate 10 f = -1420.6 |proj g|= 0.40793
## At iterate 20 f = -1420.6 |proj g|= 0.03382
##
## iterations 20
## function evaluations 25
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0338201
## final function value -1420.64
##
## F = -1420.64
## final value -1420.637053
## converged## ##
## ***************************## Estimation of CV-set 4/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1425.1 |proj g|= 17.198
## At iterate 10 f = -1425.9 |proj g|= 0.20796
##
## iterations 17
## function evaluations 30
## segments explored during Cauchy searches 17
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00811343
## final function value -1425.92
##
## F = -1425.92
## final value -1425.924348
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1161.3 |proj g|= 14
## At iterate 10 f = -1425.7 |proj g|= 0.64768
## At iterate 20 f = -1425.9 |proj g|= 0.0094272
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
##
## iterations 22
## function evaluations 55
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00522634
## final function value -1425.92
##
## F = -1425.92
## final value -1425.924355
## converged## ##
## ***************************## Estimation of CV-set 5/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1452.2 |proj g|= 9.4382
## At iterate 10 f = -1452.5 |proj g|= 0.064196
##
## iterations 18
## function evaluations 28
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0416277
## final function value -1452.46
##
## F = -1452.46
## final value -1452.455454
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1190.7 |proj g|= 14
## At iterate 10 f = -1452.4 |proj g|= 0.37119
##
## iterations 19
## function evaluations 23
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00256525
## final function value -1452.46
##
## F = -1452.46
## final value -1452.455454
## converged## ##
## ***************************## Estimation of CV-set 6/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1423.6 |proj g|= 10.14
## At iterate 10 f = -1424 |proj g|= 0.020021
##
## iterations 11
## function evaluations 29
## segments explored during Cauchy searches 12
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0200214
## final function value -1424
##
## F = -1424
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1424.001300
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1154.4 |proj g|= 14
## At iterate 10 f = -1424 |proj g|= 0.30775
##
## iterations 17
## function evaluations 34
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0714014
## final function value -1424
##
## F = -1424
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1424.001300
## converged## ##
## ***************************## Estimation of CV-set 7/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1431.9 |proj g|= 7.9183
## At iterate 10 f = -1434.1 |proj g|= 1.3391
## At iterate 20 f = -1434.3 |proj g|= 0.40425
##
## iterations 26
## function evaluations 32
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00966739
## final function value -1434.28
##
## F = -1434.28
## final value -1434.275814
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1161.7 |proj g|= 14
## At iterate 10 f = -1434.1 |proj g|= 0.40557
## At iterate 20 f = -1434.3 |proj g|= 0.3707
## At iterate 30 f = -1434.3 |proj g|= 0.072941
##
## iterations 36
## function evaluations 57
## segments explored during Cauchy searches 40
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0147916
## final function value -1434.28
##
## F = -1434.28
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1434.275785
## converged## ##
## ***************************## Estimation of CV-set 8/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1423.6 |proj g|= 7.5046
## At iterate 10 f = -1424.4 |proj g|= 2.8114
## At iterate 20 f = -1424.6 |proj g|= 0.06844
##
## iterations 26
## function evaluations 28
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00782461
## final function value -1424.61
##
## F = -1424.61
## final value -1424.606786
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1162.7 |proj g|= 14
## At iterate 10 f = -1424.1 |proj g|= 0.74494
## At iterate 20 f = -1424.6 |proj g|= 0.052911
##
## iterations 23
## function evaluations 38
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0144033
## final function value -1424.61
##
## F = -1424.61
## final value -1424.606818
## converged## ##
## ***************************## Estimation of CV-set 9/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1444.3 |proj g|= 12.803
## At iterate 10 f = -1445 |proj g|= 0.4324
## At iterate 20 f = -1445 |proj g|= 0.12333
##
## iterations 21
## function evaluations 26
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.102739
## final function value -1445.01
##
## F = -1445.01
## final value -1445.012827
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1174.4 |proj g|= 14
## At iterate 10 f = -1444.9 |proj g|= 0.4452
## At iterate 20 f = -1445 |proj g|= 0.12382
##
## iterations 22
## function evaluations 30
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0346693
## final function value -1445.01
##
## F = -1445.01
## final value -1445.012645
## converged## ##
## ***************************## Estimation of CV-set 10/10## Optimisation using starting value 1/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1449.9 |proj g|= 10.14
## At iterate 10 f = -1452.2 |proj g|= 1.4276
##
## iterations 19
## function evaluations 24
## segments explored during Cauchy searches 20
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0118335
## final function value -1452.23
##
## F = -1452.23
## final value -1452.231580
## converged## ## Optimisation using starting value 2/2## N = 5, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1180 |proj g|= 14
## At iterate 10 f = -1452.2 |proj g|= 0.75447
## At iterate 20 f = -1452.2 |proj g|= 0.0056523
##
## iterations 21
## function evaluations 26
## segments explored during Cauchy searches 25
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00106247
## final function value -1452.23
##
## F = -1452.23
## final value -1452.231549
## converged## print(est.denver.A.cv)## Cross-validation parameter estimation for STmodel
## with 10 CV-groups and 2 starting points.
## Results: 10 converged, 0 not converged.
##
## No fixed parameters.
##
## Estimated function values and convergence info:
## value convergence conv eigen.min eigen.all.min
## 1 1431.711 TRUE TRUE 0.21528208 NA
## 2 1435.925 TRUE TRUE 1.95316317 NA
## 3 1420.637 TRUE TRUE 2.37610361 NA
## 4 1425.924 TRUE TRUE 0.60717608 NA
## 5 1452.455 TRUE TRUE 1.07613640 NA
## 6 1424.001 TRUE TRUE 1.50085818 NA
## 7 1434.276 TRUE TRUE 0.06991585 NA
## 8 1424.607 TRUE TRUE 0.06886468 NA
## 9 1445.013 TRUE TRUE 0.88310030 NA
## 10 1452.232 TRUE TRUE 2.66410130 NApar(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.A, pars="all"),
plotCI((1:length(par))+.3, par, uiw=1.96*sd,
col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.A.cv, "all", boxwex=.4, col="grey", add=TRUE)
#' Save the results as an .rdata object
save(denver.data.A, denver.model.A, est.denver.model.A, est.denver.A.cv,
file = here::here("Results", "Denver_ST_Model_A.rdata"))3.1.4 Prediction using the CV model
Making predictions using the CV model. Printing out the CV summary statistics as well
pred.A.cv <- predictCV(denver.model.A, est.denver.A.cv, LTA = T)
pred.A.cv.log <- predictCV(denver.model.A, est.denver.A.cv,
LTA = T, transform="unbiased")
head(pred.A.cv$pred.all$EX)## central d_1 d_10 d_11 d_12
## 2008-12-29 NA -0.36228648 -0.232745184 -0.216254789 -0.227981974
## 2009-01-05 NA -0.28277641 -0.161511433 -0.144077302 -0.158306997
## 2009-01-12 NA -0.12950120 -0.016765871 0.002857285 -0.023815708
## 2009-01-19 NA -0.15196738 -0.046503108 -0.027625136 -0.043730230
## 2009-01-26 NA -0.09185841 0.006385409 0.025928854 0.008920688
## 2009-02-02 NA -0.01059003 0.080960068 0.101556042 0.080178246
## d_13 d_14 d_15 d_16 d_17
## 2008-12-29 -0.27587718 -0.45308211 -0.41132938 0.2175431 -0.408425796
## 2009-01-05 -0.20646305 -0.37304442 -0.32350036 0.2626777 -0.318590537
## 2009-01-12 -0.07222790 -0.22837754 -0.16069099 0.3730717 -0.159188935
## 2009-01-19 -0.09238802 -0.23853533 -0.17709189 0.3300518 -0.166740871
## 2009-01-26 -0.03996632 -0.17677830 -0.10984996 0.3611382 -0.097038051
## 2009-02-02 0.03108329 -0.09725951 -0.02162496 0.4128319 -0.008492464
## d_18 d_19 d_2 d_20 d_21
## 2008-12-29 -0.331038019 -0.29222585 0.3901875 -0.2463607208 -0.3460726863
## 2009-01-05 -0.253982498 -0.21299537 0.4284750 -0.1730029523 -0.2677362867
## 2009-01-12 -0.103754862 -0.05749019 0.5410440 -0.0249093985 -0.1160163337
## 2009-01-19 -0.127707195 -0.08346099 0.4800610 -0.0542805594 -0.1390662727
## 2009-01-26 -0.069672425 -0.02392275 0.5039757 0.0003105882 -0.0799363876
## 2009-02-02 0.009487613 0.05771830 0.5523269 0.0768787046 0.0003005563
## d_22 d_23 d_25 d_26 d_28
## 2008-12-29 -0.33307875 -0.07373230 -0.25021318 -0.27685309 -0.2974751351
## 2009-01-05 -0.25255124 -0.01185879 -0.18488305 -0.20382816 -0.2294987022
## 2009-01-12 -0.09850499 0.11497190 -0.04416628 -0.05714902 -0.0861834205
## 2009-01-19 -0.11971811 0.08771213 -0.08175036 -0.08542597 -0.1212759568
## 2009-01-26 -0.05868548 0.13350760 -0.03427905 -0.03098616 -0.0714792477
## 2009-02-02 0.02333979 0.19854577 0.03600967 0.04505800 0.0009191248
## d_29 d_3 d_31 d_32 d_33
## 2008-12-29 -0.5944721 -0.26394191 -0.03926273 -0.44919379 -0.44143570
## 2009-01-05 -0.4996518 -0.18999990 0.01415939 -0.36925205 -0.36005273
## 2009-01-12 -0.3317455 -0.05131864 0.14181773 -0.21619031 -0.20243394
## 2009-01-19 -0.3392755 -0.06721566 0.09478918 -0.23742520 -0.22637814
## 2009-01-26 -0.2656605 -0.01081514 0.13197366 -0.17685420 -0.16494842
## 2009-02-02 -0.1722824 0.06384413 0.19244453 -0.09539324 -0.08159139
## d_34 d_35 d_36 d_37 d_38
## 2008-12-29 -0.27050641 -0.4494165 -0.29216521 -0.328626549 -0.28894626
## 2009-01-05 -0.19399142 -0.3648959 -0.21965675 -0.253748051 -0.21270655
## 2009-01-12 -0.04115271 -0.2107129 -0.07389406 -0.104161228 -0.06327999
## 2009-01-19 -0.06968022 -0.2232687 -0.10212755 -0.132100589 -0.08800042
## 2009-01-26 -0.01252819 -0.1582361 -0.04808845 -0.076173118 -0.03068253
## 2009-02-02 0.06694804 -0.0739274 0.02744664 0.001607336 0.04782714
## d_39 d_4 d_40 d_41 d_42
## 2008-12-29 -0.34337745 -0.30736245 -0.5345659 -0.36635499 -0.368793538
## 2009-01-05 -0.27222325 -0.22171269 -0.4364337 -0.29292362 -0.283314659
## 2009-01-12 -0.12715302 -0.02968994 -0.2321539 -0.14475780 -0.091459697
## 2009-01-19 -0.15748644 -0.08872076 -0.2794322 -0.17405966 -0.150651402
## 2009-01-26 -0.10472010 -0.02746285 -0.2072056 -0.11940384 -0.089543652
## 2009-02-02 -0.03011308 0.06648047 -0.1033111 -0.04277704 0.004263438
## d_43 d_44 d_45 d_46 d_47
## 2008-12-29 -0.35374440 -0.39385807 -0.305318718 -0.28919264 -0.27786268
## 2009-01-05 -0.27927994 -0.31137824 -0.235347298 -0.21503560 -0.20388790
## 2009-01-12 -0.13159658 -0.15541494 -0.090578943 -0.06565131 -0.05354352
## 2009-01-19 -0.15798845 -0.17478990 -0.123138425 -0.09492467 -0.08446268
## 2009-01-26 -0.10223054 -0.11204170 -0.071522991 -0.03969683 -0.02954283
## 2009-02-02 -0.02513611 -0.02846003 0.002345503 0.03762877 0.04790834
## d_48 d_49 d_5 d_50 d_51
## 2008-12-29 -0.28784879 -0.241872686 -0.373821107 -0.326813616 -0.5147255
## 2009-01-05 -0.21420855 -0.172601938 -0.289902488 -0.253699278 -0.4279910
## 2009-01-12 -0.06533173 -0.038507580 -0.132298341 -0.105338871 -0.2716342
## 2009-01-19 -0.09509167 -0.058802693 -0.150613781 -0.135593958 -0.2821055
## 2009-01-26 -0.04031796 -0.006506988 -0.086630858 -0.081282376 -0.2151275
## 2009-02-02 0.03659565 0.064428316 -0.001847949 -0.004788027 -0.1290538
## d_52 d_53 d_54 d_55 d_56
## 2008-12-29 -0.28711260 -0.12261843 -0.27544123 -0.5808377 -0.366916130
## 2009-01-05 -0.21192210 -0.05280246 -0.20218427 -0.4841291 -0.284728369
## 2009-01-12 -0.06202891 0.12367243 -0.05254475 -0.2812472 -0.128823825
## 2009-01-19 -0.08967451 0.04973375 -0.08413976 -0.3298659 -0.148768905
## 2009-01-26 -0.03347287 0.09707791 -0.02985068 -0.2588902 -0.086306954
## 2009-02-02 0.04455631 0.17839841 0.04702823 -0.1561306 -0.002903898
## d_57 d_58 d_6 d_60 d_61
## 2008-12-29 -0.5430766 -0.29262610 -0.323992416874 -0.30174711 -0.31409306
## 2009-01-05 -0.4525655 -0.21521568 -0.250514279586 -0.22655601 -0.23298237
## 2009-01-12 -0.2888908 -0.06440499 -0.103390117553 -0.07501728 -0.07853822
## 2009-01-19 -0.3004780 -0.08832676 -0.131240370363 -0.10479126 -0.09897610
## 2009-01-26 -0.2306498 -0.03001057 -0.076402303282 -0.04880258 -0.03740834
## 2009-02-02 -0.1407071 0.04948817 0.000003149964 0.02961824 0.04504029
## d_62 d_63 d_64 d_65 d_66
## 2008-12-29 -0.41729878 -0.3975971 -0.5041282 -0.39328711 -0.27193223
## 2009-01-05 -0.32935511 -0.3157525 -0.4205455 -0.30840497 -0.20124275
## 2009-01-12 -0.16802663 -0.1641974 -0.2624116 -0.15049204 -0.05662886
## 2009-01-19 -0.18225728 -0.1792727 -0.2821558 -0.16707543 -0.08739982
## 2009-01-26 -0.11470780 -0.1165916 -0.2185796 -0.10216311 -0.03504185
## 2009-02-02 -0.02677031 -0.0344163 -0.1338600 -0.01676362 0.03919469
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59
## 2008-12-29 -0.19999832 -0.35203089 -0.27310791 -0.30107370 NA NA NA NA
## 2009-01-05 -0.12616073 -0.26923058 -0.19730288 -0.22400540 NA NA NA NA
## 2009-01-12 0.05426315 -0.11673689 -0.04806857 -0.07714022 NA NA NA NA
## 2009-01-19 -0.01588908 -0.13091238 -0.07350184 -0.09671253 NA NA NA NA
## 2009-01-26 0.03498902 -0.06739144 -0.01659636 -0.03822853 NA NA NA NA
## 2009-02-02 0.11951558 0.01554583 0.06162256 0.04013914 NA NA NA NA
## d_9
## 2008-12-29 NA
## 2009-01-05 NA
## 2009-01-12 NA
## 2009-01-19 NA
## 2009-01-26 NA
## 2009-02-02 NAtail(pred.A.cv$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 1.0153075 1.0365885 1.0403129 1.0482544 1.0144463 0.9451461
## 2020-11-09 NA 0.6129213 0.6313776 0.6211391 0.6343057 0.5830190 0.5403579
## 2020-11-16 NA 0.9616024 0.9988670 0.9817190 1.0003479 0.9577870 0.8917466
## 2020-11-23 NA 0.8771148 0.9064736 0.8988584 0.9086319 0.8676386 0.7912111
## 2020-11-30 NA 1.0243965 1.0570188 1.0569792 1.0663244 1.0334303 0.9434250
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2
## 2020-11-02 1.1085251 1.1489975 1.1710023 1.0253529 1.0865084 1.2657413
## 2020-11-09 0.7044141 0.7257865 0.7560962 0.6175884 0.6846931 0.8208796
## 2020-11-16 1.0735201 1.1230581 1.1124280 0.9765002 1.0537491 1.2079120
## 2020-11-23 0.9723052 1.0133818 1.0188022 0.8865214 0.9578493 1.1278449
## 2020-11-30 1.1125346 1.2048538 1.1681545 1.0376216 1.1021347 1.2913014
## 2020-12-07 NA NA NA NA NA NA
## d_20 d_21 d_22 d_23 d_25 d_26
## 2020-11-02 1.0200958 1.0374512 1.0716942 1.0775568 0.9168185 1.0357455
## 2020-11-09 0.6135664 0.6255490 0.6636695 0.6716347 0.5135114 0.6045204
## 2020-11-16 0.9891965 0.9871666 1.0176236 1.0523952 0.8828499 0.9770651
## 2020-11-23 0.8986404 0.8977975 0.9308337 0.9539043 0.7897345 0.8918454
## 2020-11-30 1.0689908 1.0486209 1.0800445 1.1080207 0.9462193 1.0531385
## 2020-12-07 NA NA NA NA NA NA
## d_28 d_29 d_3 d_31 d_32 d_33
## 2020-11-02 0.8933947 1.0406468 1.0939037 0.9419153 0.9546268 0.9500936
## 2020-11-09 0.4985135 0.6348037 0.6717188 0.5426915 0.5477225 0.5550487
## 2020-11-16 0.8561805 0.9841973 1.0426941 0.9181546 0.9064500 0.8995412
## 2020-11-23 0.7661558 0.8892249 0.9487463 0.8287602 0.8137409 0.8143944
## 2020-11-30 0.9178555 1.0339336 1.1074392 0.9826598 0.9642513 0.9577187
## 2020-12-07 NA NA NA NA NA NA
## d_34 d_35 d_36 d_37 d_38 d_39
## 2020-11-02 1.0762004 1.0494578 0.9920703 0.9937540 1.0610605 0.8902132
## 2020-11-09 0.6642520 0.6338048 0.5864955 0.5850722 0.6497705 0.5031871
## 2020-11-16 1.0304952 0.9968617 0.9476558 0.9502138 1.0150518 0.8487846
## 2020-11-23 0.9387189 0.8999141 0.8576758 0.8581069 0.9195882 0.7640492
## 2020-11-30 1.0948636 1.0555747 1.0097651 1.0102803 1.0762677 0.9086520
## 2020-12-07 NA NA NA NA NA NA
## d_4 d_40 d_41 d_42 d_43 d_44
## 2020-11-02 1.0395736 0.8535904 0.9384733 0.9905828 0.9672294 1.0466225
## 2020-11-09 0.6278772 0.4747605 0.5277814 0.5831185 0.5605568 0.6346899
## 2020-11-16 0.9808899 0.7888230 0.8972960 0.9403258 0.9247074 0.9905022
## 2020-11-23 0.9094250 0.6893282 0.7999222 0.8633025 0.8293803 0.9028538
## 2020-11-30 1.0531131 0.8155815 0.9586894 1.0047726 0.9844266 1.0529555
## 2020-12-07 NA NA NA NA NA NA
## d_45 d_46 d_47 d_48 d_49 d_5
## 2020-11-02 0.9303928 0.9916843 1.0129575 0.9969385 1.0382311 1.0273759
## 2020-11-09 0.5281338 0.5936733 0.6169624 0.5962082 0.6248963 0.6316690
## 2020-11-16 0.8909502 0.9422248 0.9830546 0.9516216 0.9999965 0.9672184
## 2020-11-23 0.7991530 0.8602363 0.8929020 0.8655292 0.9037291 0.8789137
## 2020-11-30 0.9525905 1.0073778 1.0289339 1.0139954 1.0571023 1.0225070
## 2020-12-07 NA NA NA NA NA NA
## d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 0.9280294 1.0091265 1.0366260 1.1157954 0.9993623 0.9486490
## 2020-11-09 0.5403765 0.6048557 0.6228188 0.6716793 0.5925683 0.5497244
## 2020-11-16 0.8832730 0.9679981 0.9849256 1.0494501 0.9637032 0.8978995
## 2020-11-23 0.8012276 0.8704512 0.8982857 0.9753801 0.8749914 0.8151933
## 2020-11-30 0.9445205 1.0101989 1.0516888 1.1248030 1.0389529 0.9508642
## 2020-12-07 NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61
## 2020-11-02 1.0291476 1.0143207 1.0992731 0.9797674 0.9921893 1.0747313
## 2020-11-09 0.6269554 0.6146639 0.6820153 0.5707252 0.5870810 0.6667059
## 2020-11-16 0.9685909 0.9605498 1.0482364 0.9351738 0.9428642 1.0134782
## 2020-11-23 0.8847239 0.8682101 0.9582192 0.8446043 0.8522978 0.9295177
## 2020-11-30 1.0316212 1.0118100 1.1078981 0.9969283 1.0040087 1.0795107
## 2020-12-07 NA NA NA NA NA NA
## d_62 d_63 d_64 d_65 d_66 d_67
## 2020-11-02 1.0705972 1.0484585 0.9326397 1.0942474 0.9797545 0.9938950
## 2020-11-09 0.6681986 0.6388058 0.5362505 0.6543527 0.5715386 0.5896571
## 2020-11-16 1.0063436 0.9954318 0.8829836 1.0208160 0.9361592 0.9660635
## 2020-11-23 0.9218063 0.9041185 0.7933088 0.9373491 0.8440889 0.8821866
## 2020-11-30 1.0666864 1.0544009 0.9377339 1.1003707 1.0005123 1.0210878
## 2020-12-07 NA NA NA NA NA NA
## d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.1566223 1.0348270 1.0718938 NA NA NA NA NA
## 2020-11-09 0.7310982 0.6245033 0.6607406 NA NA NA NA NA
## 2020-11-16 1.0958777 0.9774879 1.0209648 NA NA NA NA NA
## 2020-11-23 1.0038295 0.8910882 0.9309178 NA NA NA NA NA
## 2020-11-30 1.1551023 1.0442159 1.0832408 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA NA NAhead(pred.A.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 0.7142497 0.8123136 0.8269625 0.8179457 0.7796934 0.6530775
## 2009-01-05 NA 0.7732058 0.8722277 0.8886006 0.8766751 0.8354583 0.7072605
## 2009-01-12 NA 0.9011275 1.0080113 1.0289849 1.0025813 0.9552003 0.8171081
## 2009-01-19 NA 0.8809779 0.9784260 0.9978760 0.9825665 0.9359015 0.8086473
## 2009-01-26 NA 0.9354427 1.0315220 1.0525857 1.0354722 0.9860685 0.8599840
## 2009-02-02 NA 1.0145411 1.1113505 1.1351168 1.1117679 1.0585039 0.9310091
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2008-12-29 0.6793596 1.277067 0.6818200 0.7365967 0.7656016 1.518451 0.8024370
## 2009-01-05 0.7417020 1.335580 0.7457778 0.7954396 0.8286148 1.577126 0.8632654
## 2009-01-12 0.8728193 1.491027 0.8745214 0.9242156 0.9679108 1.764458 1.0008064
## 2009-01-19 0.8586003 1.427886 0.8678293 0.9022054 0.9430001 1.659605 0.9716286
## 2009-01-26 0.9183014 1.472667 0.9303781 0.9559953 1.0007652 1.699380 1.0259628
## 2009-02-02 1.0029841 1.550543 1.0164302 1.0346460 1.0858237 1.783240 1.1074480
## d_21 d_22 d_23 d_25 d_26 d_28
## 2008-12-29 0.7252808 0.7366993 0.9543645 0.7981263 0.7793073 0.7612828
## 2009-01-05 0.7843267 0.7981790 1.0149399 0.8519815 0.8380318 0.8148053
## 2009-01-12 0.9127702 0.9308048 1.1518457 0.9806869 0.9701062 0.9403349
## 2009-01-19 0.8919249 0.9110118 1.1205897 0.9444900 0.9427942 0.9078866
## 2009-01-26 0.9462141 0.9681217 1.1728591 0.9903878 0.9953129 0.9542220
## 2009-02-02 1.0252295 1.0506875 1.2514702 1.0624890 1.0737561 1.0258522
## d_29 d_3 d_31 d_32 d_33 d_34
## 2008-12-29 0.5657541 0.7890550 0.9865918 0.6545086 0.6594803 0.7824119
## 2009-01-05 0.6219814 0.8493264 1.0405254 0.7088368 0.7152971 0.8445123
## 2009-01-12 0.7356504 0.9753831 1.1820029 0.8259299 0.8373115 0.9838539
## 2009-01-19 0.7300936 0.9597593 1.1275350 0.8084543 0.8174169 0.9560853
## 2009-01-26 0.7858334 1.0152366 1.1701087 0.8588302 0.8691315 1.0122337
## 2009-02-02 0.8627188 1.0937557 1.2429293 0.9316282 0.9446216 1.0958920
## d_35 d_36 d_37 d_38 d_39 d_4
## 2008-12-29 0.6544368 0.7657941 0.7390663 0.7682631 0.7278840 0.7547770
## 2009-01-05 0.7120317 0.8232176 0.7963009 0.8289591 0.7814088 0.8220648
## 2009-01-12 0.8306037 0.9522296 0.9245522 0.9623904 0.9032459 0.9958726
## 2009-01-19 0.8201335 0.9255812 0.8988836 0.9387498 0.8761291 0.9386090
## 2009-01-26 0.8751480 0.9768538 0.9504191 0.9940057 0.9234883 0.9977483
## 2009-02-02 0.9520495 1.0533951 1.0271489 1.0750841 0.9949263 1.0958884
## d_40 d_41 d_42 d_43 d_44 d_45
## 2008-12-29 0.6013752 0.7117019 0.7098057 0.7200597 0.6932568 0.7564946
## 2009-01-05 0.6632144 0.7657085 0.7729522 0.7755694 0.7525790 0.8110891
## 2009-01-12 0.8133450 0.8877702 0.9362192 0.8988387 0.8793119 0.9371954
## 2009-01-19 0.7756392 0.8619476 0.8822437 0.8752950 0.8621972 0.9069758
## 2009-01-26 0.8336039 0.9102072 0.9376908 0.9253712 0.9178203 0.9548490
## 2009-02-02 0.9247552 0.9825564 1.0297833 0.9994358 0.9976476 1.0279073
## d_46 d_47 d_48 d_49 d_5 d_50
## 2008-12-29 0.7676143 0.7766774 0.7686466 0.8066624 0.7060585 0.7392724
## 2009-01-05 0.8266755 0.8361958 0.8273595 0.8642322 0.7677155 0.7953232
## 2009-01-12 0.9598415 0.9717383 0.9601483 0.9879592 0.8986104 0.9224938
## 2009-01-19 0.9321286 0.9420559 0.9319729 0.9678678 0.8821712 0.8949801
## 2009-01-26 0.9850366 0.9951566 0.9844249 1.0196199 0.9403456 0.9449133
## 2009-02-02 1.0642106 1.0752239 1.0631117 1.0943948 1.0234492 1.0200142
## d_51 d_52 d_53 d_54 d_55 d_56
## 2008-12-29 0.6130619 0.7703936 0.9079292 0.7785604 0.5741824 0.7109506
## 2009-01-05 0.6684939 0.8303133 0.9733361 0.8376216 0.6323246 0.7716980
## 2009-01-12 0.7815128 0.9643380 1.1609358 0.9727094 0.7743795 0.9017381
## 2009-01-19 0.7732716 0.9378403 1.0779900 0.9423602 0.7374910 0.8838002
## 2009-01-26 0.8267494 0.9918811 1.1300779 0.9948502 0.7916139 0.9406502
## 2009-02-02 0.9009868 1.0722249 1.2256647 1.0742780 0.8771777 1.0223690
## d_57 d_58 d_6 d_60 d_61 d_62
## 2008-12-29 0.5955915 0.7650992 0.7434237 0.7583467 0.7488499 0.6771953
## 2009-01-05 0.6519686 0.8266209 0.7998065 0.8174541 0.8120643 0.7391708
## 2009-01-12 0.7678618 0.9611164 0.9262688 0.9510939 0.9476283 0.8682919
## 2009-01-19 0.7589761 0.9383486 0.9005751 0.9230986 0.9284088 0.8557829
## 2009-01-26 0.8138332 0.9946538 0.9511207 0.9761735 0.9873227 0.9153765
## 2009-02-02 0.8903940 1.0769189 1.0264518 1.0557367 1.0721395 0.9993348
## d_63 d_64 d_65 d_66 d_67 d_68
## 2008-12-29 0.6892434 0.6201032 0.6921422 0.7817906 0.8403231 0.7213761
## 2009-01-05 0.7478973 0.6739660 0.7533067 0.8388895 0.9044899 0.7835130
## 2009-01-12 0.8701524 0.7892311 0.8820142 0.9692466 1.0830890 0.9124459
## 2009-01-19 0.8570216 0.7736336 0.8673772 0.9397370 1.0095204 0.8994860
## 2009-01-26 0.9123626 0.8242687 0.9254336 0.9901302 1.0620464 0.9583736
## 2009-02-02 0.9904190 0.8970170 1.0078387 1.0663282 1.1555778 1.0411594
## d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.7801793 0.7590881 NA NA NA NA NA
## 2009-01-05 0.8415614 0.8197610 NA NA NA NA NA
## 2009-01-12 0.9769465 0.9493006 NA NA NA NA NA
## 2009-01-19 0.9523631 0.9307804 NA NA NA NA NA
## 2009-01-26 1.0080862 0.9867341 NA NA NA NA NA
## 2009-02-02 1.0900662 1.0670824 NA NA NA NA NAtail(pred.A.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 2.781967 2.834168 2.850092 2.871261 2.773636 2.590173
## 2020-11-09 NA 1.860271 1.889887 1.874058 1.897848 1.801553 1.727810
## 2020-11-16 NA 2.636266 2.729150 2.687548 2.736549 2.620466 2.455126
## 2020-11-23 NA 2.422608 2.488260 2.473722 2.496600 2.394446 2.220186
## 2020-11-30 NA 2.806957 2.892492 2.897385 2.922914 2.826112 2.585098
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2020-11-02 3.042630 3.168553 3.244899 2.804743 2.977606 3.565357 2.787609
## 2020-11-09 2.031147 2.075041 2.142842 1.865439 1.992255 2.284875 1.856295
## 2020-11-16 2.937909 3.086961 3.060044 2.670773 2.881450 3.364473 2.702442
## 2020-11-23 2.655084 2.766155 2.786469 2.440876 2.617895 3.105413 2.468360
## 2020-11-30 3.054758 3.349759 3.235245 2.838943 3.024171 3.656686 2.926678
## 2020-12-07 NA NA NA NA NA NA NA
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2020-11-02 2.838464 2.944943 2.950946 2.514945 2.836488 2.459586 2.849619
## 2020-11-09 1.880129 1.958100 1.966235 1.680234 1.842735 1.657154 1.898990
## 2020-11-16 2.699164 2.789486 2.877194 2.430905 2.674396 2.369693 2.693118
## 2020-11-23 2.468375 2.557450 2.607189 2.214754 2.455786 2.165668 2.449085
## 2020-11-30 2.870175 2.968854 3.041489 2.589904 2.885488 2.520415 2.830385
## 2020-12-07 NA NA NA NA NA NA NA
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2020-11-02 3.003200 2.578411 2.612859 2.606174 2.949157 2.871160 2.712563
## 2020-11-09 1.968774 1.729614 1.739312 1.755585 1.953326 1.894619 1.808085
## 2020-11-16 2.852857 2.517638 2.489736 2.477540 2.817212 2.723829 2.594486
## 2020-11-23 2.596908 2.302268 2.269219 2.275265 2.570108 2.472086 2.371153
## 2020-11-30 3.043392 2.685236 2.637735 2.625844 3.004389 2.888393 2.760578
## 2020-12-07 NA NA NA NA NA NA NA
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2020-11-02 2.717994 2.903920 2.455849 2.851992 2.365285 2.570651 2.713591
## 2020-11-09 1.806047 1.924605 1.667618 1.889399 1.619319 1.704710 1.805334
## 2020-11-16 2.601857 2.773088 2.355974 2.689141 2.216697 2.466630 2.580294
## 2020-11-23 2.372806 2.520521 2.164496 2.503570 2.006685 2.237667 2.388918
## 2020-11-30 2.762705 2.947977 2.501186 2.890340 2.276651 2.622596 2.751865
## 2020-12-07 NA NA NA NA NA NA NA
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2020-11-02 2.644784 2.871290 2.552176 2.716652 2.767308 2.728546 2.838609
## 2020-11-09 1.760972 1.901682 1.706792 1.824634 1.862357 1.827645 1.877416
## 2020-11-16 2.534449 2.714152 2.453156 2.585510 2.685604 2.607608 2.731720
## 2020-11-23 2.303931 2.486247 2.237889 2.381970 2.454025 2.392490 2.480876
## 2020-11-30 2.690259 2.888769 2.608915 2.759540 2.811569 2.775402 2.891986
## 2020-12-07 NA NA NA NA NA NA NA
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 2.816501 2.550457 2.756640 2.838729 3.068558 2.730138 2.603335
## 2020-11-09 1.895985 1.730845 1.839873 1.876630 1.968020 1.817608 1.746836
## 2020-11-16 2.651826 2.438780 2.645348 2.695350 2.871252 2.634324 2.474236
## 2020-11-23 2.427622 2.246664 2.399420 2.471544 2.666157 2.410638 2.277745
## 2020-11-30 2.802412 2.592789 2.759231 2.881209 3.095746 2.840086 2.608630
## 2020-12-07 NA NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 2.821899 2.777065 3.017327 2.683165 2.713423 2.950928 2.945523
## 2020-11-09 1.887339 1.862124 1.987928 1.782228 1.809528 1.962216 1.969537
## 2020-11-16 2.655848 2.631589 2.867090 2.565725 2.582660 2.775499 2.761772
## 2020-11-23 2.442118 2.399440 2.620244 2.343427 2.358982 2.551949 2.537752
## 2020-11-30 2.828475 2.769939 3.043280 2.728884 2.745389 2.964892 2.933258
## 2020-12-07 NA NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2020-11-02 2.871315 2.561443 3.004188 2.679337 2.715772 3.195663 2.833098
## 2020-11-09 1.906125 1.723075 1.934912 1.781227 1.812622 2.088038 1.879540
## 2020-11-16 2.722806 2.437048 2.791237 2.564807 2.640928 3.007077 2.675123
## 2020-11-23 2.485121 2.227918 2.567636 2.339134 2.428353 2.742560 2.453666
## 2020-11-30 2.888048 2.573989 3.022190 2.735129 2.790114 3.190385 2.859654
## 2020-12-07 NA NA NA NA NA NA NA
## d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.938774 NA NA NA NA NA
## 2020-11-09 1.947983 NA NA NA NA NA
## 2020-11-16 2.792628 NA NA NA NA NA
## 2020-11-23 2.552077 NA NA NA NA NA
## 2020-11-30 2.971919 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NAsummary(pred.A.cv)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX
## obs 0.16756831 0.16756831 0.10939771
## average 0.09377047 0.09377047 0.05212872
##
## R2:
## EX.mu EX.mu.beta EX
## obs 0.5831563 0.5831563 0.8223333
## average 0.6093873 0.6093873 0.8792830
##
## Coverage of 95% prediction intervals:
## EX
## obs 0.9443038
## average 0.8571429summary(pred.A.cv.log)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.2182824 0.2182824 0.15097844 0.15098707
## average 0.1045726 0.1045726 0.06003835 0.06050165
##
## R2:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.6108816 0.6108816 0.8138452 0.8138239
## average 0.6750348 0.6750348 0.8928829 0.8912233
##
## Coverage of 95% prediction intervals:
## EX.pred
## obs 0.9443038
## average 0.9047619par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModA.jpeg"),
width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.A.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.A.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()## quartz_off_screen
## 23.1.5 Predictions for 2009-2020
What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites
x.A <- coef(est.denver.model.A, pars = "cov")$par
x.A## [1] -7.797037 -7.526946 12.449157 -3.008551 -4.860112E.A <- predict(denver.model.A, est.denver.model.A, LTA = T,
transform="unbiased", pred.var=FALSE)
print(E.A)## Prediction for STmodel.
##
## Regression parameters:
## 0 Spatio-temporal covariate(s).
## 28 beta-fields regression parameters in x$pars.
##
## Regression parameters are assumed to be known and
## prediction variances do NOT include
## uncertainties in regression parameters.
##
## Prediction of beta-fields, (x$beta):
## List of 2
## $ mu: num [1:69, 1:2] 0.2557 0.028 0.0441 0.0712 0.042 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX: num [1:69, 1:2] 0.2646 0.0409 0.0405 0.062 0.0302 ...
## ..- attr(*, "dimnames")=List of 2
##
## Predictions for log-Gaussian field of type: unbiased
##
## Predictions for 622 times at 69 locations.
## List of 5
## $ EX.mu : num [1:622, 1:69] 1.21 1.27 1.43 1.36 1.41 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.mu.beta: num [1:622, 1:69] 1.26 1.32 1.48 1.41 1.45 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX : num [1:622, 1:69] 1.29 1.35 1.52 1.44 1.48 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.pred : num [1:622, 1:69] 1.3 1.36 1.53 1.45 1.49 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.EX : num [1:622, 1:69] 0.231 0.276 0.394 0.341 0.37 ...
## ..- attr(*, "dimnames")=List of 2
##
## Variances have been computed.
## List of 2
## $ log.VX : num [1:622, 1:69] 0.0499 0.0498 0.0497 0.0497 0.0496 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.VX.pred: num [1:622, 1:69] 0.0576 0.0575 0.0575 0.0574 0.0574 ...
## ..- attr(*, "dimnames")=List of 2
##
## Mean squared prediciton errors NOT been computed.
##
## 69 temporal averages have been compute.
## List of 1
## $ LTA:'data.frame': 69 obs. of 4 variables:week_preds <- as.data.frame(E.A$EX) %>%
mutate(week = as.Date(rownames(E.A$EX)),
year = year(as.Date(rownames(E.A$EX))))
week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]
box_preds <- week_preds %>%
select(week, all_of(week_sites)) %>%
#filter(week %in% week_weeks) %>%
mutate(month = month(week),
year = year(week)) %>%
filter(year %in% c(2009:2020))
box_data <- box_preds %>%
pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)## week month year location
## Min. :2009-01-05 Min. : 1.000 Min. :2009 Length:42228
## 1st Qu.:2012-01-09 1st Qu.: 4.000 1st Qu.:2012 Class :character
## Median :2014-12-29 Median : 7.000 Median :2014 Mode :character
## Mean :2014-12-28 Mean : 6.494 Mean :2014
## 3rd Qu.:2017-12-18 3rd Qu.: 9.000 3rd Qu.:2017
## Max. :2020-12-07 Max. :12.000 Max. :2020
##
## pred
## Min. :0.3202
## 1st Qu.:1.1202
## Median :1.3738
## Mean :1.4175
## 3rd Qu.:1.6484
## Max. :3.8293
## NA's :861#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
group_by(month) %>%
summarize(median = median(pred, na.rm=T)) %>%
arrange(desc(median))
box_summarypred_box_plot <- ggplot(box_data) +
geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
show.legend = F) +
#scale_color_viridis(name = "Month", discrete = T) +
facet_wrap(. ~ year, scales = "free") +
xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
# theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot
3.2 Model B: 2 temporal trend (basis) functions
3.2.1 Create the model object
For this version of the model, use iid for cov.beta (beta, beta1, beta2) and exp for cov.nu (error). Here we can specify different LUR formluae. The length of the LUR list should be number of basis functions + 1.
denver.data.B <- denver.data2.2
names(denver.data.B$covars)## [1] "ID" "lon" "lat" "tree_cover_100"
## [5] "impervious_2500" "low_int_100" "med_int_50" "med_int_2500"
## [9] "high_int_100" "pop_den_50" "dist_m_cafos" "dist_m_og_wells"
## [13] "dist_m_wwtp" "aadt_100" "aadt_250" "aadt_1000"
## [17] "x" "y" "type"LUR.B <- list(covar_fun, covar_fun, covar_fun)
cov.beta.B <- list(covf = c("iid", "iid", "iid"), nugget = T)
cov.nu.B <- list(covf = "exp", nugget = T, random.effect = FALSE)
locations.B <- list(coords = c("x","y"), long.lat = c("lon","lat"), others= "type")
denver.model.B <- createSTmodel(denver.data.B, LUR = LUR.B,
ST = c("bc_st_no2", "bc_st_smk"),
cov.beta = cov.beta.B, cov.nu = cov.nu.B,
locations = locations.B)## No trend $trend.fnc object detected, STdata probably from old version of the package.
## $trend.fnc has been added based on spline fit to elements in STmodel$trend.denver.model.B## STmodel-object with:
## No. locations: 69 (observed: 64)
## No. time points: 622 (observed: 231)
## No. obs: 1021
##
## Trend with 2 basis function(s):
## [1] "V1" "V2"
## with dates:
## 2008-12-29 to 2020-12-07
##
## Models for the beta-fields are:
## $const
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos +
## dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## $V1
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos +
## dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## $V2
## ~~tree_cover_100 + impervious_2500 + low_int_100 + med_int_50 +
## med_int_2500 + high_int_100 + pop_den_50 + dist_m_cafos +
## dist_m_og_wells + dist_m_wwtp + aadt_100 + aadt_250 + aadt_1000
##
## 2 spatio-temporal covariate(s):
## [1] "bc_st_no2" "bc_st_smk"
##
## Covariance model for the beta-field(s):
## Covariance type(s): iid, iid, iid
## Nugget: Yes, Yes, Yes
## Covariance model for the nu-field(s):
## Covariance type: exp
## Nugget: ~1
## Random effect: No
## All sites:
## central dist
## 1 68
## Observed:
## central dist
## 1 63
##
## For central:
## Number of obs: 231
## Dates: 2016-02-08 to 2020-11-30
## For dist:
## Number of obs: 790
## Dates: 2018-05-07 to 2020-04-063.2.2 Estimate model parameters
- nugget values can range from -1 to -6
- range values can range from 0 to 4
- try starting values where the sill and nugget are the same (e.g., both 0) and where they are very different (e.g., 0 and -5)
- make sure the starting values are pretty different
names <- loglikeSTnames(denver.model.B, all=FALSE)
names## [1] "log.nugget.const.iid" "log.nugget.V1.iid"
## [3] "log.nugget.V2.iid" "nu.log.range.exp"
## [5] "nu.log.sill.exp" "nu.log.nugget.(Intercept).exp"# x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
# c(-1, -1, -1, 0, -1, -1),
# c(-1, -1, -1, 0, -5, -1),
# c(-5, -5, -5, 0, -1, -5),
# c(-5, -5, -5, 0, -5, -5),
# c(-1, -1, -1, 2, -1, -1),
# c(-1, -1, -1, 2, -5, -1),
# c(-5, -5, -5, 2, -1, -5),
# c(-5, -5, -5, 2, -5, -5),
# c(-1, -1, -1, 4, -1, -1),
# c(-1, -1, -1, 4, -5, -1),
# c(-5, -5, -5, 4, -1, -5),
# c(-5, -5, -5, 4, -5, -5))
x.init.B <- cbind(c(0, 0, 0, 0, 0, 0),
c(-10, -5, -5, 4, -3, -5),
c(-10, -5, -5, 4, -5, -5),
c(-10, -5, -5, 6, -3, -5),
c(-10, -5, -5, 6, -5, -5),
c(-10, -5, -5, 8, -3, -5),
c(-10, -5, -5, 8, -5, -5),
c(-10, -5, -5, 10, -3, -5),
c(-10, -5, -5, 10, -5, -5),
c(-12, -5, -5, 8, -5, -5),
c(-12, -5, -5, 10, -5, -5),
c(-12, -5, -5, 12, -5, -5))
rownames(x.init.B) <- loglikeSTnames(denver.model.B, all=FALSE)
x.init.B## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## log.nugget.const.iid 0 -10 -10 -10 -10 -10 -10 -10 -10
## log.nugget.V1.iid 0 -5 -5 -5 -5 -5 -5 -5 -5
## log.nugget.V2.iid 0 -5 -5 -5 -5 -5 -5 -5 -5
## nu.log.range.exp 0 4 4 6 6 8 8 10 10
## nu.log.sill.exp 0 -3 -5 -3 -5 -3 -5 -3 -5
## nu.log.nugget.(Intercept).exp 0 -5 -5 -5 -5 -5 -5 -5 -5
## [,10] [,11] [,12]
## log.nugget.const.iid -12 -12 -12
## log.nugget.V1.iid -5 -5 -5
## log.nugget.V2.iid -5 -5 -5
## nu.log.range.exp 8 10 12
## nu.log.sill.exp -5 -5 -5
## nu.log.nugget.(Intercept).exp -5 -5 -5#' Difference from tutorial: use Josh Keller's version of the function
source(here::here("Code", "functions_model.R"))
est.denver.model.B <- estimate.STmodel(denver.model.B, x.init.B)## Optimisation using starting value 1/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= 507.58 |proj g|= 15
## At iterate 10 f = -1393.9 |proj g|= 2.5132
## At iterate 20 f = -1394.2 |proj g|= 0.013016
##
## iterations 21
## function evaluations 35
## segments explored during Cauchy searches 26
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0130169
## final function value -1394.21
##
## F = -1394.21
## final value -1394.206501
## converged## Optimisation using starting value 2/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1197.7 |proj g|= 12
## At iterate 10 f = -1393.8 |proj g|= 3.6232
## At iterate 20 f = -1394.1 |proj g|= 2.0253
## ys=-8.565e+00 -gs= 2.159e+00, BFGS update SKIPPED
## At iterate 30 f = -1608.5 |proj g|= 20.131
##
## iterations 38
## function evaluations 58
## segments explored during Cauchy searches 41
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.168749
## final function value -1622.14
##
## F = -1622.14
## final value -1622.137478
## converged## Optimisation using starting value 3/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1237.6 |proj g|= 20
## At iterate 10 f = -1392.3 |proj g|= 3.6831
## At iterate 20 f = -1394.1 |proj g|= 0.28437
## ys=-5.332e-02 -gs= 1.204e+00, BFGS update SKIPPED
## At iterate 30 f = -1606.3 |proj g|= 18.399
## At iterate 40 f = -1622.1 |proj g|= 0.40807
##
## iterations 45
## function evaluations 68
## segments explored during Cauchy searches 49
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0219034
## final function value -1622.13
##
## F = -1622.13
## final value -1622.130756
## converged## Optimisation using starting value 4/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1202.8 |proj g|= 12
## At iterate 10 f = -1631.2 |proj g|= 7.6237
## At iterate 20 f = -1635.5 |proj g|= 0.79559
## At iterate 30 f = -1635.5 |proj g|= 0.11189
## At iterate 40 f = -1635.5 |proj g|= 0.090055
##
## iterations 43
## function evaluations 62
## segments explored during Cauchy searches 47
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00445772
## final function value -1635.49
##
## F = -1635.49
## final value -1635.494460
## converged## Optimisation using starting value 5/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1245.1 |proj g|= 20
## At iterate 10 f = -1613.5 |proj g|= 3.2561
## At iterate 20 f = -1624.6 |proj g|= 7.3182
## At iterate 30 f = -1626.9 |proj g|= 4.5415
## At iterate 40 f = -1630.8 |proj g|= 1.2802
## At iterate 50 f = -1631.1 |proj g|= 1.002
## ys=-4.238e-01 -gs= 3.760e-01, BFGS update SKIPPED
## At iterate 60 f = -1635.3 |proj g|= 1.4564
## At iterate 70 f = -1635.5 |proj g|= 0.35894
## At iterate 80 f = -1635.5 |proj g|= 0.2711
##
## iterations 88
## function evaluations 136
## segments explored during Cauchy searches 93
## BFGS updates skipped 1
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0246402
## final function value -1635.49
##
## F = -1635.49
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1635.494327
## converged## Optimisation using starting value 6/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1262.8 |proj g|= 12
## At iterate 10 f = -1635.3 |proj g|= 3.9677
## At iterate 20 f = -1635.5 |proj g|= 0.072351
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1635.492358
## stopped after 28 iterations## Optimisation using starting value 7/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1297.5 |proj g|= 20
## At iterate 10 f = -1634.5 |proj g|= 2.4806
## At iterate 20 f = -1635.5 |proj g|= 0.011806
## At iterate 30 f = -1635.5 |proj g|= 0.063142
## At iterate 40 f = -1635.5 |proj g|= 0.0048559
##
## iterations 41
## function evaluations 69
## segments explored during Cauchy searches 46
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00485598
## final function value -1635.49
##
## F = -1635.49
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1635.494445
## converged## Optimisation using starting value 8/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1459.9 |proj g|= 12
## At iterate 10 f = -1635.4 |proj g|= 2.6314
## At iterate 20 f = -1635.5 |proj g|= 0.029661
##
## iterations 22
## function evaluations 26
## segments explored during Cauchy searches 27
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00905988
## final function value -1635.49
##
## F = -1635.49
## final value -1635.492637
## converged## Optimisation using starting value 9/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1370.4 |proj g|= 20
## At iterate 10 f = -1635.5 |proj g|= 0.85406
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1635.492327
## stopped after 14 iterations## Optimisation using starting value 10/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1297.3 |proj g|= 20
## At iterate 10 f = -1634.9 |proj g|= 2.363
##
## iterations 19
## function evaluations 42
## segments explored during Cauchy searches 24
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0598606
## final function value -1635.49
##
## F = -1635.49
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1635.491845
## converged## Optimisation using starting value 11/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1370.2 |proj g|= 20
## At iterate 10 f = -1635.4 |proj g|= 1.0976
##
## iterations 15
## function evaluations 21
## segments explored during Cauchy searches 20
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00186049
## final function value -1635.49
##
## F = -1635.49
## final value -1635.491833
## converged## Optimisation using starting value 12/12## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1373.5 |proj g|= 20
## At iterate 10 f = -1635.5 |proj g|= 0.19953
##
## iterations 15
## function evaluations 31
## segments explored during Cauchy searches 19
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0461308
## final function value -1635.49
##
## F = -1635.49
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1635.491840
## convergedprint(est.denver.model.B)## Optimisation for STmodel with 12 starting points.
## Results: 4 converged, 8 not converged, 0 failed.
## Best result for starting point 4, optimisation has converged
##
## No fixed parameters.
##
## Estimated parameters for all starting point(s):
## [,1] [,2] [,3]
## gamma.bc_st_no2 0.0060118836 0.00596005607 0.005955238637
## gamma.bc_st_smk 0.0521824404 0.05195395316 0.051963119374
## alpha.const.(Intercept) 0.1912785980 0.25013512306 0.250232500893
## alpha.const.tree_cover_100 0.0032725799 0.00858563403 0.008600000285
## alpha.const.impervious_2500 0.0086346267 0.01913402640 0.019148742535
## alpha.const.low_int_100 0.0093760168 0.01430191351 0.014301194760
## alpha.const.med_int_50 -0.0039613701 -0.00421505227 -0.004210495232
## alpha.const.med_int_2500 0.0171963057 -0.00094776366 -0.000943197470
## alpha.const.high_int_100 0.0014554826 0.01229113752 0.012293428681
## alpha.const.pop_den_50 -0.0158271837 -0.01315055513 -0.013153454350
## alpha.const.dist_m_cafos -0.0196381194 -0.02064697856 -0.020650426238
## alpha.const.dist_m_og_wells -0.0014823607 0.00410607533 0.004094448529
## alpha.const.dist_m_wwtp 0.0061514019 0.00280361913 0.002820074602
## alpha.const.aadt_100 0.0344456466 0.01913393117 0.019130992908
## alpha.const.aadt_250 0.0216541261 0.02327931343 0.023287934150
## alpha.const.aadt_1000 -0.0082169471 -0.00884375176 -0.008843479197
## alpha.V1.(Intercept) 0.1640948824 0.18448242880 0.184486002539
## alpha.V1.tree_cover_100 -0.0149371422 -0.01216815574 -0.012157621514
## alpha.V1.impervious_2500 0.0046700767 -0.00589806167 -0.005876364566
## alpha.V1.low_int_100 -0.0099590774 -0.00751872533 -0.007518516324
## alpha.V1.med_int_50 -0.0145381299 -0.01258187668 -0.012580731420
## alpha.V1.med_int_2500 0.0007319780 0.00362165651 0.003607057482
## alpha.V1.high_int_100 0.0193619626 0.01805976475 0.018074702267
## alpha.V1.pop_den_50 0.0027890733 0.00271971642 0.002727439623
## alpha.V1.dist_m_cafos 0.0059170892 -0.00001866864 -0.000001210275
## alpha.V1.dist_m_og_wells -0.0016883212 0.00379434591 0.003797105516
## alpha.V1.dist_m_wwtp -0.0027537848 -0.00835900122 -0.008348554976
## alpha.V1.aadt_100 -0.0244963327 -0.03129890664 -0.031308202681
## alpha.V1.aadt_250 -0.0206482343 -0.00905272933 -0.009062198058
## alpha.V1.aadt_1000 0.0044713009 0.00044568561 0.000442089992
## alpha.V2.(Intercept) 0.1604857126 0.18506394002 0.185100939006
## alpha.V2.tree_cover_100 -0.0038598984 -0.00140760762 -0.001396523996
## alpha.V2.impervious_2500 -0.0292271566 -0.02688644737 -0.026870279271
## alpha.V2.low_int_100 -0.0008278762 0.00081574734 0.000813085668
## alpha.V2.med_int_50 -0.0085606701 -0.00852572347 -0.008520145449
## alpha.V2.med_int_2500 0.0213995519 0.01609592640 0.016090534357
## alpha.V2.high_int_100 -0.0115300275 -0.00531046255 -0.005307474586
## alpha.V2.pop_den_50 0.0070030441 0.00627161540 0.006274912103
## alpha.V2.dist_m_cafos 0.0143984442 0.01517681980 0.015169579556
## alpha.V2.dist_m_og_wells -0.0018465363 0.00316452234 0.003156172950
## alpha.V2.dist_m_wwtp 0.0009020415 -0.00309514261 -0.003081324416
## alpha.V2.aadt_100 -0.0189470731 -0.02803138572 -0.028037993207
## alpha.V2.aadt_250 0.0107539488 0.01242409203 0.012422670722
## alpha.V2.aadt_1000 -0.0056846619 -0.00812879018 -0.008124952189
## log.nugget.const.iid -14.2075145954 -14.99853953428 -14.998066086823
## log.nugget.V1.iid -13.3885208673 -13.25977435576 -13.811919503385
## log.nugget.V2.iid -15.0000000000 -14.99612179953 -14.990498086304
## nu.log.range.exp 0.0000000000 13.75456682509 13.762410630096
## nu.log.sill.exp -4.6993271432 -3.00667112679 -3.006919056540
## nu.log.nugget.(Intercept).exp -4.2038352592 -4.76493602901 -4.764763069367
## [,4] [,5] [,6]
## gamma.bc_st_no2 0.0062756184 0.0062750888 0.0062448184
## gamma.bc_st_smk 0.0518645345 0.0518633379 0.0518689417
## alpha.const.(Intercept) 0.2490535778 0.2490551735 0.2494632803
## alpha.const.tree_cover_100 0.0048752209 0.0048744305 0.0049646999
## alpha.const.impervious_2500 0.0196052934 0.0196284557 0.0197915602
## alpha.const.low_int_100 0.0139807586 0.0139840762 0.0140120667
## alpha.const.med_int_50 -0.0090675828 -0.0090833998 -0.0091188953
## alpha.const.med_int_2500 -0.0042252843 -0.0042602214 -0.0044315319
## alpha.const.high_int_100 0.0098436326 0.0098418272 0.0098587736
## alpha.const.pop_den_50 -0.0097277268 -0.0097230856 -0.0097633578
## alpha.const.dist_m_cafos -0.0212607744 -0.0212768491 -0.0213592746
## alpha.const.dist_m_og_wells 0.0063053930 0.0063061812 0.0062501046
## alpha.const.dist_m_wwtp -0.0014448112 -0.0014503791 -0.0014192174
## alpha.const.aadt_100 0.0198837181 0.0198758207 0.0198041276
## alpha.const.aadt_250 0.0232382337 0.0232435263 0.0233246181
## alpha.const.aadt_1000 -0.0102055575 -0.0102042173 -0.0101736236
## alpha.V1.(Intercept) 0.1796396617 0.1796034928 0.1797166847
## alpha.V1.tree_cover_100 -0.0095807532 -0.0095857651 -0.0095930706
## alpha.V1.impervious_2500 -0.0026502125 -0.0026663502 -0.0026861356
## alpha.V1.low_int_100 -0.0063794936 -0.0063868015 -0.0064147407
## alpha.V1.med_int_50 -0.0115235700 -0.0115134344 -0.0114282183
## alpha.V1.med_int_2500 -0.0001957307 -0.0001779400 -0.0000751741
## alpha.V1.high_int_100 0.0210467647 0.0210293676 0.0209677001
## alpha.V1.pop_den_50 0.0030558366 0.0030481096 0.0029986132
## alpha.V1.dist_m_cafos -0.0014346281 -0.0014369914 -0.0014559278
## alpha.V1.dist_m_og_wells 0.0011285833 0.0011302192 0.0012182994
## alpha.V1.dist_m_wwtp -0.0045033269 -0.0045077458 -0.0045269901
## alpha.V1.aadt_100 -0.0254282735 -0.0254231814 -0.0254880303
## alpha.V1.aadt_250 -0.0107957045 -0.0107872255 -0.0107660633
## alpha.V1.aadt_1000 -0.0003513547 -0.0003443569 -0.0003134788
## alpha.V2.(Intercept) 0.1837226965 0.1837249799 0.1838253301
## alpha.V2.tree_cover_100 -0.0032709767 -0.0032725068 -0.0032224787
## alpha.V2.impervious_2500 -0.0224524331 -0.0224434950 -0.0223627890
## alpha.V2.low_int_100 0.0020186127 0.0020206665 0.0020292072
## alpha.V2.med_int_50 -0.0113583401 -0.0113640954 -0.0113460420
## alpha.V2.med_int_2500 0.0096917376 0.0096847193 0.0096683086
## alpha.V2.high_int_100 -0.0053929727 -0.0053929997 -0.0053821449
## alpha.V2.pop_den_50 0.0087220406 0.0087151215 0.0086414050
## alpha.V2.dist_m_cafos 0.0125478935 0.0125360189 0.0124589374
## alpha.V2.dist_m_og_wells 0.0035149242 0.0035251114 0.0035775837
## alpha.V2.dist_m_wwtp -0.0032472277 -0.0032520674 -0.0032526984
## alpha.V2.aadt_100 -0.0273235923 -0.0273252702 -0.0273759384
## alpha.V2.aadt_250 0.0125332866 0.0125343813 0.0125506316
## alpha.V2.aadt_1000 -0.0082224575 -0.0082276268 -0.0082511959
## log.nugget.const.iid -10.4632900888 -10.3985631570 -9.9979407159
## log.nugget.V1.iid -7.0285457395 -7.0243931348 -7.0258773232
## log.nugget.V2.iid -7.7309663192 -7.7317201739 -7.7474076123
## nu.log.range.exp 13.5125634739 13.5098031295 13.5109546812
## nu.log.sill.exp -2.9784990884 -2.9783780489 -2.9780490360
## nu.log.nugget.(Intercept).exp -4.9340326598 -4.9344286674 -4.9352276245
## [,7] [,8] [,9]
## gamma.bc_st_no2 0.0062750305 0.00624754768 0.00624517200
## gamma.bc_st_smk 0.0518629046 0.05186815364 0.05186799148
## alpha.const.(Intercept) 0.2490496581 0.24942294174 0.24945251190
## alpha.const.tree_cover_100 0.0048788656 0.00495761154 0.00496455640
## alpha.const.impervious_2500 0.0196042830 0.01977807254 0.01979078798
## alpha.const.low_int_100 0.0139810141 0.01400994065 0.01401226517
## alpha.const.med_int_50 -0.0090654055 -0.00911525896 -0.00911852466
## alpha.const.med_int_2500 -0.0042249847 -0.00441729279 -0.00443215946
## alpha.const.high_int_100 0.0098442364 0.00985798842 0.00985905250
## alpha.const.pop_den_50 -0.0097324016 -0.00976044755 -0.00976428759
## alpha.const.dist_m_cafos -0.0212606928 -0.02135274314 -0.02135967695
## alpha.const.dist_m_og_wells 0.0063029166 0.00625448630 0.00625008615
## alpha.const.dist_m_wwtp -0.0014443857 -0.00142087749 -0.00141947406
## alpha.const.aadt_100 0.0198820425 0.01981001382 0.01980396759
## alpha.const.aadt_250 0.0232404112 0.02331764231 0.02332424882
## alpha.const.aadt_1000 -0.0102036016 -0.01017580725 -0.01017292717
## alpha.V1.(Intercept) 0.1796525981 0.17970823595 0.17971635560
## alpha.V1.tree_cover_100 -0.0095817610 -0.00959386626 -0.00959521360
## alpha.V1.impervious_2500 -0.0026517806 -0.00268701723 -0.00269063387
## alpha.V1.low_int_100 -0.0063784085 -0.00641321983 -0.00641533893
## alpha.V1.med_int_50 -0.0115203868 -0.01143558341 -0.01142777810
## alpha.V1.med_int_2500 -0.0001905905 -0.00008132414 -0.00007091681
## alpha.V1.high_int_100 0.0210459615 0.02097097055 0.02096485396
## alpha.V1.pop_den_50 0.0030534950 0.00300235265 0.00299744584
## alpha.V1.dist_m_cafos -0.0014387133 -0.00145376537 -0.00145682140
## alpha.V1.dist_m_og_wells 0.0011348363 0.00121167969 0.00121976868
## alpha.V1.dist_m_wwtp -0.0045044776 -0.00452732568 -0.00452936877
## alpha.V1.aadt_100 -0.0254339827 -0.02548402101 -0.02548962933
## alpha.V1.aadt_250 -0.0107952680 -0.01076688982 -0.01076441074
## alpha.V1.aadt_1000 -0.0003501826 -0.00031577454 -0.00031259535
## alpha.V2.(Intercept) 0.1837201618 0.18381566669 0.18382257167
## alpha.V2.tree_cover_100 -0.0032693538 -0.00322654205 -0.00322298325
## alpha.V2.impervious_2500 -0.0224532625 -0.02237127152 -0.02236530503
## alpha.V2.low_int_100 0.0020185719 0.00202819139 0.00202903497
## alpha.V2.med_int_50 -0.0113557183 -0.01134716779 -0.01134570891
## alpha.V2.med_int_2500 0.0096941148 0.00967168608 0.00967055008
## alpha.V2.high_int_100 -0.0053933123 -0.00538316208 -0.00538260579
## alpha.V2.pop_den_50 0.0087181759 0.00864682945 0.00863992996
## alpha.V2.dist_m_cafos 0.0125471507 0.01246616230 0.01245955972
## alpha.V2.dist_m_og_wells 0.0035164838 0.00357346862 0.00357858028
## alpha.V2.dist_m_wwtp -0.0032481932 -0.00325304734 -0.00325403892
## alpha.V2.aadt_100 -0.0273253755 -0.02737162877 -0.02737577550
## alpha.V2.aadt_250 0.0125341258 0.01254927405 0.01255077775
## alpha.V2.aadt_1000 -0.0082225089 -0.00824950170 -0.00825156424
## log.nugget.const.iid -10.4481949545 -10.02637857089 -9.99617297360
## log.nugget.V1.iid -7.0300419363 -7.02605082295 -7.02622038251
## log.nugget.V2.iid -7.7315355639 -7.74665482316 -7.74794236783
## nu.log.range.exp 13.5121556541 13.51065516608 13.51030811048
## nu.log.sill.exp -2.9784997305 -2.97815642597 -2.97807591297
## nu.log.nugget.(Intercept).exp -4.9340481552 -4.93516126191 -4.93527648027
## [,10] [,11] [,12]
## gamma.bc_st_no2 0.0063171711 0.0063168583 0.0063113077
## gamma.bc_st_smk 0.0518590513 0.0518631206 0.0518749984
## alpha.const.(Intercept) 0.2485173757 0.2485413475 0.2486812126
## alpha.const.tree_cover_100 0.0047559128 0.0047561336 0.0047672588
## alpha.const.impervious_2500 0.0193399493 0.0193531093 0.0193839872
## alpha.const.low_int_100 0.0139381060 0.0139393168 0.0139435416
## alpha.const.med_int_50 -0.0089886588 -0.0089936373 -0.0089963190
## alpha.const.med_int_2500 -0.0039321095 -0.0039414894 -0.0039561026
## alpha.const.high_int_100 0.0098228588 0.0098238091 0.0098299831
## alpha.const.pop_den_50 -0.0096881146 -0.0096810069 -0.0096742703
## alpha.const.dist_m_cafos -0.0211182900 -0.0211236440 -0.0211321894
## alpha.const.dist_m_og_wells 0.0063769747 0.0063774050 0.0063714233
## alpha.const.dist_m_wwtp -0.0014812408 -0.0014771607 -0.0014596871
## alpha.const.aadt_100 0.0199933751 0.0199919591 0.0199874581
## alpha.const.aadt_250 0.0231216219 0.0231228210 0.0231301077
## alpha.const.aadt_1000 -0.0102441085 -0.0102462022 -0.0102482989
## alpha.V1.(Intercept) 0.1795708992 0.1795570796 0.1795671190
## alpha.V1.tree_cover_100 -0.0095625102 -0.0095602429 -0.0095541761
## alpha.V1.impervious_2500 -0.0025956207 -0.0025908271 -0.0025747943
## alpha.V1.low_int_100 -0.0063244373 -0.0063289545 -0.0063375917
## alpha.V1.med_int_50 -0.0116509561 -0.0116532175 -0.0116568871
## alpha.V1.med_int_2500 -0.0003597753 -0.0003655688 -0.0003772143
## alpha.V1.high_int_100 0.0211614029 0.0211601862 0.0211612946
## alpha.V1.pop_den_50 0.0031321987 0.0031348258 0.0031406720
## alpha.V1.dist_m_cafos -0.0014137351 -0.0014040404 -0.0013821293
## alpha.V1.dist_m_og_wells 0.0010156107 0.0010083986 0.0010010728
## alpha.V1.dist_m_wwtp -0.0044677063 -0.0044667258 -0.0044660557
## alpha.V1.aadt_100 -0.0253585159 -0.0253526404 -0.0253508028
## alpha.V1.aadt_250 -0.0108364819 -0.0108368589 -0.0108408111
## alpha.V1.aadt_1000 -0.0004058246 -0.0004074424 -0.0004118690
## alpha.V2.(Intercept) 0.1835930005 0.1836037312 0.1836445815
## alpha.V2.tree_cover_100 -0.0033363676 -0.0033349391 -0.0033253831
## alpha.V2.impervious_2500 -0.0225745467 -0.0225672019 -0.0225518680
## alpha.V2.low_int_100 0.0020053069 0.0020054011 0.0020050406
## alpha.V2.med_int_50 -0.0113689933 -0.0113722966 -0.0113732523
## alpha.V2.med_int_2500 0.0097242073 0.0097193770 0.0097146054
## alpha.V2.high_int_100 -0.0054088191 -0.0054068723 -0.0054014877
## alpha.V2.pop_den_50 0.0088287099 0.0088324299 0.0088371225
## alpha.V2.dist_m_cafos 0.0126710514 0.0126685543 0.0126643029
## alpha.V2.dist_m_og_wells 0.0034284932 0.0034273196 0.0034236107
## alpha.V2.dist_m_wwtp -0.0032402629 -0.0032366303 -0.0032262324
## alpha.V2.aadt_100 -0.0272558331 -0.0272558188 -0.0272598777
## alpha.V2.aadt_250 0.0125113744 0.0125105768 0.0125097559
## alpha.V2.aadt_1000 -0.0081800766 -0.0081813330 -0.0081832037
## log.nugget.const.iid -11.9828653970 -11.9963559520 -12.0001162114
## log.nugget.V1.iid -7.0355605905 -7.0331561647 -7.0306406707
## log.nugget.V2.iid -7.7081729398 -7.7085438681 -7.7106093212
## nu.log.range.exp 13.5151755391 13.5167035286 13.5224235635
## nu.log.sill.exp -2.9785879051 -2.9790274847 -2.9789261257
## nu.log.nugget.(Intercept).exp -4.9320850646 -4.9321253281 -4.9314759991
##
## Function value(s):
## [1] 1394.207 1622.137 1622.131 1635.494 1635.494 1635.492 1635.494 1635.493
## [9] 1635.492 1635.492 1635.492 1635.4923.2.3 Cross-validation
Define the CV groups
set.seed(1000)
site_idsB <- colnames(denver.model.B$D.beta)[which(str_detect(colnames(denver.model.B$D.beta), "d_") == T)]
site_idsB## [1] "d_1" "d_10" "d_11" "d_12" "d_13" "d_14" "d_15" "d_16" "d_17" "d_18"
## [11] "d_19" "d_2" "d_20" "d_21" "d_22" "d_23" "d_25" "d_26" "d_28" "d_29"
## [21] "d_3" "d_31" "d_32" "d_33" "d_34" "d_35" "d_36" "d_37" "d_38" "d_39"
## [31] "d_4" "d_40" "d_41" "d_42" "d_43" "d_44" "d_45" "d_46" "d_47" "d_48"
## [41] "d_49" "d_5" "d_50" "d_51" "d_52" "d_53" "d_54" "d_55" "d_56" "d_57"
## [51] "d_58" "d_6" "d_60" "d_61" "d_62" "d_63" "d_64" "d_65" "d_66" "d_67"
## [61] "d_68" "d_7" "d_8"Ind.cv.B <- createCV(denver.model.B, groups = 10, #min.dist = .1,
subset = site_idsB)
ID.cv.B <- sapply(split(denver.model.B$obs$ID, Ind.cv.B), unique)
print(sapply(ID.cv.B, length))## 0 1 2 3 4 5 6 7 8 9 10
## 1 7 7 7 6 6 6 6 6 6 6table(Ind.cv.B)## Ind.cv.B
## 0 1 2 3 4 5 6 7 8 9 10
## 231 75 80 97 73 64 93 80 81 66 81I.col.B <- apply(sapply(ID.cv.B, function(x) denver.model.B$locations$ID%in% x), 1,
function(x) if(sum(x)==1) which(x) else 0)
names(I.col.B) <- denver.model.B$locations$ID
print(I.col.B)## central d_1 d_10 d_11 d_12 d_13 d_14 d_15 d_16 d_17
## 1 6 2 3 4 4 4 9 4 10
## d_18 d_19 d_2 d_20 d_21 d_22 d_23 d_25 d_26 d_28
## 8 5 11 3 2 11 4 9 11 9
## d_29 d_3 d_31 d_32 d_33 d_34 d_35 d_36 d_37 d_38
## 2 4 6 8 5 5 10 8 3 8
## d_39 d_4 d_40 d_41 d_42 d_43 d_44 d_45 d_46 d_47
## 6 7 7 3 7 8 11 3 9 5
## d_48 d_49 d_5 d_50 d_51 d_52 d_53 d_54 d_55 d_56
## 9 4 6 9 10 3 7 5 7 6
## d_57 d_58 d_6 d_60 d_61 d_62 d_63 d_64 d_65 d_66
## 2 2 11 5 2 11 10 3 8 6
## d_67 d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 7 10 2 10 0 0 0 0 0par(mfrow=c(1,1))
plot(denver.model.B$locations$long,
denver.model.B$locations$lat,
pch=23+floor(I.col.B/max(I.col.B)+.5), bg=I.col.B,
xlab="Longitude", ylab="Latitude")
map("county", "colorado", col="#FFFF0055",fill=TRUE, add=TRUE)## [[1]]
## NULLID starting conditions using previous model without CV:
x.init.B.cv <- coef(est.denver.model.B, pars="cov")[,c("par","init")]
x.init.B.cvRun the model with cross validation
est.denver.B.cv <- estimateCV(denver.model.B, x.init.B.cv, Ind.cv.B)##
## ***************************## Estimation of CV-set 1/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1484.4 |proj g|= 13.035
## At iterate 10 f = -1485.3 |proj g|= 0.076883
## At iterate 20 f = -1485.3 |proj g|= 0.0093172
##
## iterations 21
## function evaluations 35
## segments explored during Cauchy searches 21
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00931736
## final function value -1485.34
##
## F = -1485.34
## final value -1485.339955
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1105.4 |proj g|= 12
## At iterate 10 f = -1479.4 |proj g|= 0.73963
## At iterate 20 f = -1483.4 |proj g|= 13.853
## At iterate 30 f = -1485.3 |proj g|= 0.94886
##
## iterations 39
## function evaluations 64
## segments explored during Cauchy searches 42
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00705453
## final function value -1485.34
##
## F = -1485.34
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1485.343747
## converged## ##
## ***************************## Estimation of CV-set 2/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1500 |proj g|= 8.8707
## At iterate 10 f = -1500.3 |proj g|= 0.39016
## At iterate 20 f = -1500.3 |proj g|= 0.24659
## At iterate 30 f = -1500.4 |proj g|= 0.5606
## At iterate 40 f = -1500.4 |proj g|= 0.071363
##
## iterations 45
## function evaluations 54
## segments explored during Cauchy searches 46
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00240616
## final function value -1500.38
##
## F = -1500.38
## final value -1500.383951
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1104.7 |proj g|= 12
## At iterate 10 f = -1424.8 |proj g|= 11.786
## At iterate 20 f = -1487.1 |proj g|= 0.93927
## At iterate 30 f = -1489.7 |proj g|= 2.6261
## At iterate 40 f = -1490.3 |proj g|= 1.4331
## At iterate 50 f = -1494 |proj g|= 1.6876
## At iterate 60 f = -1494.5 |proj g|= 0.13324
## At iterate 70 f = -1494.5 |proj g|= 0.10735
## At iterate 80 f = -1499.1 |proj g|= 4.7371
## At iterate 90 f = -1500.3 |proj g|= 0.79519
## At iterate 100 f = -1500.4 |proj g|= 0.0523
##
## iterations 104
## function evaluations 137
## segments explored during Cauchy searches 107
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.143867
## final function value -1500.38
##
## F = -1500.38
## final value -1500.382816
## converged## ##
## ***************************## Estimation of CV-set 3/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1467.3 |proj g|= 10.066
## At iterate 10 f = -1467.6 |proj g|= 0.15397
##
## iterations 14
## function evaluations 18
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0369096
## final function value -1467.62
##
## F = -1467.62
## final value -1467.616601
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1089.1 |proj g|= 12
## At iterate 10 f = -1400 |proj g|= 12.375
## At iterate 20 f = -1458 |proj g|= 2.2193
## At iterate 30 f = -1460.4 |proj g|= 0.14393
## At iterate 40 f = -1463.9 |proj g|= 1.0706
## At iterate 50 f = -1463.9 |proj g|= 0.013697
## At iterate 60 f = -1466.9 |proj g|= 1.7069
## At iterate 70 f = -1467.6 |proj g|= 0.45824
## At iterate 80 f = -1467.6 |proj g|= 0.1933
##
## iterations 82
## function evaluations 122
## segments explored during Cauchy searches 85
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0399222
## final function value -1467.62
##
## F = -1467.62
## final value -1467.618744
## converged## ##
## ***************************## Estimation of CV-set 4/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1481.7 |proj g|= 17.346
## At iterate 10 f = -1482.6 |proj g|= 0.3585
## At iterate 20 f = -1482.6 |proj g|= 0.032367
## At iterate 30 f = -1482.7 |proj g|= 0.36671
## At iterate 40 f = -1482.7 |proj g|= 0.12535
## At iterate 50 f = -1482.7 |proj g|= 0.39405
##
## iterations 56
## function evaluations 78
## segments explored during Cauchy searches 56
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0209832
## final function value -1482.67
##
## F = -1482.67
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1482.672168
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1107.5 |proj g|= 12
## At iterate 10 f = -1416.2 |proj g|= 11.607
## At iterate 20 f = -1473.2 |proj g|= 2.2818
## At iterate 30 f = -1475.1 |proj g|= 9.4907
## At iterate 40 f = -1475.7 |proj g|= 0.02604
## At iterate 50 f = -1477.7 |proj g|= 1.0248
## At iterate 60 f = -1477.8 |proj g|= 0.07492
## At iterate 70 f = -1482.2 |proj g|= 2.4183
## At iterate 80 f = -1482.5 |proj g|= 0.80487
## At iterate 90 f = -1482.6 |proj g|= 0.37974
## At iterate 100 f = -1482.7 |proj g|= 0.38035
## At iterate 110 f = -1482.7 |proj g|= 0.057702
##
## iterations 119
## function evaluations 168
## segments explored during Cauchy searches 122
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0145089
## final function value -1482.67
##
## F = -1482.67
## final value -1482.672380
## converged## ##
## ***************************## Estimation of CV-set 5/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1506 |proj g|= 10.656
## At iterate 10 f = -1506.5 |proj g|= 0.059847
## At iterate 20 f = -1506.6 |proj g|= 0.3413
## At iterate 30 f = -1506.6 |proj g|= 0.15811
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1506.581892
## stopped after 35 iterations## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1120.3 |proj g|= 12
## At iterate 10 f = -1493.7 |proj g|= 12.155
## At iterate 20 f = -1499.8 |proj g|= 4.429
## At iterate 30 f = -1506.3 |proj g|= 1.8316
## At iterate 40 f = -1506.6 |proj g|= 0.18166
## At iterate 50 f = -1506.6 |proj g|= 0.39657
## At iterate 60 f = -1506.6 |proj g|= 0.098565
##
## iterations 65
## function evaluations 97
## segments explored during Cauchy searches 68
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.102392
## final function value -1506.58
##
## F = -1506.58
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1506.582174
## converged## ##
## ***************************## Estimation of CV-set 6/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1481.1 |proj g|= 10.066
## At iterate 10 f = -1481.5 |proj g|= 0.06454
## At iterate 20 f = -1481.6 |proj g|= 0.90862
## At iterate 30 f = -1481.6 |proj g|= 0.07344
## Bad direction in the line search;
## refresh the lbfgs memory and restart the iteration.
## Line search cannot locate an adequate point after 20 function
## and gradient evaluations
## final value -1481.582531
## stopped after 35 iterations## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1086.4 |proj g|= 12
## At iterate 10 f = -1471.7 |proj g|= 3.8336
## At iterate 20 f = -1479.3 |proj g|= 2.7325
## At iterate 30 f = -1481.3 |proj g|= 5.9125
## At iterate 40 f = -1481.5 |proj g|= 0.33309
## At iterate 50 f = -1481.6 |proj g|= 0.25228
## At iterate 60 f = -1481.6 |proj g|= 0.034288
## At iterate 70 f = -1481.6 |proj g|= 0.027797
##
## iterations 74
## function evaluations 90
## segments explored during Cauchy searches 77
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0336023
## final function value -1481.58
##
## F = -1481.58
## final value -1481.582484
## converged## ##
## ***************************## Estimation of CV-set 7/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1491.1 |proj g|= 2.1306
## At iterate 10 f = -1492.4 |proj g|= 0.95479
## At iterate 20 f = -1492.4 |proj g|= 0.12235
## At iterate 30 f = -1492.5 |proj g|= 0.46536
## At iterate 40 f = -1492.5 |proj g|= 0.33101
## At iterate 50 f = -1492.5 |proj g|= 0.0067843
##
## iterations 50
## function evaluations 57
## segments explored during Cauchy searches 51
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.00678434
## final function value -1492.55
##
## F = -1492.55
## final value -1492.549537
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1099.7 |proj g|= 12
## At iterate 10 f = -1484.8 |proj g|= 1.04
## At iterate 20 f = -1486.8 |proj g|= 7.7877
## At iterate 30 f = -1489.2 |proj g|= 0.57899
## At iterate 40 f = -1492.4 |proj g|= 0.33942
## At iterate 50 f = -1492.5 |proj g|= 0.3765
##
## iterations 59
## function evaluations 90
## segments explored during Cauchy searches 62
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0306232
## final function value -1492.55
##
## F = -1492.55
## final value -1492.549474
## converged## ##
## ***************************## Estimation of CV-set 8/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1479.5 |proj g|= 9.3568
## At iterate 10 f = -1480.6 |proj g|= 1.2495
## At iterate 20 f = -1480.7 |proj g|= 0.48288
## At iterate 30 f = -1481 |proj g|= 0.188
## At iterate 40 f = -1481 |proj g|= 0.0058638
##
## iterations 41
## function evaluations 61
## segments explored during Cauchy searches 41
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.00586381
## final function value -1481
##
## F = -1481
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1481.002261
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1097.9 |proj g|= 12
## At iterate 10 f = -1473.3 |proj g|= 12.094
## At iterate 20 f = -1478.8 |proj g|= 2.3195
## At iterate 30 f = -1480 |proj g|= 0.67451
## At iterate 40 f = -1481 |proj g|= 0.15932
##
## iterations 49
## function evaluations 77
## segments explored during Cauchy searches 52
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.016169
## final function value -1481
##
## F = -1481
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1481.002269
## converged## ##
## ***************************## Estimation of CV-set 9/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1500.8 |proj g|= 11.701
## At iterate 10 f = -1501.3 |proj g|= 0.13133
##
## iterations 15
## function evaluations 31
## segments explored during Cauchy searches 15
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.10318
## final function value -1501.33
##
## F = -1501.33
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1501.329327
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1117.3 |proj g|= 12
## At iterate 10 f = -1493.4 |proj g|= 6.7681
## At iterate 20 f = -1497.9 |proj g|= 14.57
## At iterate 30 f = -1501.2 |proj g|= 0.35627
## At iterate 40 f = -1501.3 |proj g|= 0.62611
## At iterate 50 f = -1501.3 |proj g|= 0.036213
##
## iterations 56
## function evaluations 80
## segments explored during Cauchy searches 59
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0143426
## final function value -1501.33
##
## F = -1501.33
## final value -1501.334466
## converged## ##
## ***************************## Estimation of CV-set 10/10## Optimisation using starting value 1/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1499.1 |proj g|= 10.066
## At iterate 10 f = -1500.3 |proj g|= 0.2129
## At iterate 20 f = -1500.3 |proj g|= 0.059763
## At iterate 30 f = -1500.3 |proj g|= 0.20161
## At iterate 40 f = -1500.3 |proj g|= 0.058964
## At iterate 50 f = -1500.3 |proj g|= 0.00997
##
## iterations 52
## function evaluations 68
## segments explored during Cauchy searches 53
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 1
## norm of the final projected gradient 0.0152949
## final function value -1500.3
##
## F = -1500.3
## Warning: more than 10 function and gradient evaluations
## in the last line search
## final value -1500.301602
## converged## ## Optimisation using starting value 2/2## N = 6, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= -1099.5 |proj g|= 12
## At iterate 10 f = -1488.4 |proj g|= 5.2544
## At iterate 20 f = -1491.7 |proj g|= 5.0281
## At iterate 30 f = -1493.6 |proj g|= 10.355
## At iterate 40 f = -1498 |proj g|= 1.1003
## At iterate 50 f = -1500 |proj g|= 0.78914
## At iterate 60 f = -1500.3 |proj g|= 0.9166
## At iterate 70 f = -1500.3 |proj g|= 0.31545
## At iterate 80 f = -1500.3 |proj g|= 0.043618
##
## iterations 83
## function evaluations 128
## segments explored during Cauchy searches 86
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0248369
## final function value -1500.3
##
## F = -1500.3
## final value -1500.301353
## converged## print(est.denver.B.cv)## Cross-validation parameter estimation for STmodel
## with 10 CV-groups and 2 starting points.
## Results: 7 converged, 3 not converged.
##
## No fixed parameters.
##
## Estimated function values and convergence info:
## value convergence conv eigen.min eigen.all.min
## 1 1485.344 TRUE TRUE 0.0036525174 NA
## 2 1500.384 TRUE FALSE -0.0022060915 NA
## 3 1467.619 TRUE TRUE 0.0028793325 NA
## 4 1482.672 TRUE FALSE -0.0001855608 NA
## 5 1506.582 TRUE TRUE 0.0002412768 NA
## 6 1481.582 TRUE TRUE 0.0001981437 NA
## 7 1492.550 TRUE TRUE 0.0009942168 NA
## 8 1481.002 TRUE TRUE 0.5570127079 NA
## 9 1501.334 TRUE TRUE 0.0439724903 NA
## 10 1500.302 TRUE FALSE -0.0008568626 NApar(mfrow=c(1,1), mar=c(13.5,2.5,.5,.5), las=2)
with(coef(est.denver.model.B, pars="all"),
plotCI((1:length(par))+.3, par, uiw=1.96*sd,
col=2, xlab="", xaxt="n", ylab=""))
boxplot(est.denver.B.cv, "all", boxwex=.4, col="grey", add=TRUE)
#' Save the results as an .rdata object
save(denver.data.B, denver.model.B, est.denver.model.B, est.denver.B.cv,
file = here::here("Results", "Denver_ST_Model_B.rdata"))3.2.4 Prediction using the CV model
Making predictions using the CV model. Printing out the CV summary statistics as well
pred.B.cv <- predictCV(denver.model.B, est.denver.B.cv, LTA = T)
pred.B.cv.log <- predictCV(denver.model.B, est.denver.B.cv,
LTA = T, transform="unbiased")
head(pred.B.cv$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 0.2221147 0.3406085 0.4326163 0.4042291 0.3058418 0.3641883
## 2009-01-05 NA 0.2362779 0.3500038 0.4331282 0.4063972 0.3093446 0.3606838
## 2009-01-12 NA 0.2901571 0.3954783 0.4669599 0.4364097 0.3407006 0.3848774
## 2009-01-19 NA 0.2460560 0.3517986 0.4190743 0.3967714 0.3024217 0.3391215
## 2009-01-26 NA 0.2449216 0.3470935 0.4075001 0.3868296 0.2938626 0.3226113
## 2009-02-02 NA 0.2523968 0.3500981 0.4036145 0.3821287 0.2905731 0.3107640
## d_15 d_16 d_17 d_18 d_19 d_2
## 2008-12-29 0.1678977 0.4593548 0.2136028 0.3918066 0.2466937 0.3625573
## 2009-01-05 0.1907173 0.4715380 0.2343124 0.3924895 0.2662171 0.3934257
## 2009-01-12 0.2479878 0.5120979 0.2890843 0.4273435 0.3281190 0.4650666
## 2009-01-19 0.2179013 0.4841166 0.2578214 0.3777103 0.2854634 0.4396387
## 2009-01-26 0.2239366 0.4875188 0.2624378 0.3649943 0.2886366 0.4599268
## 2009-02-02 0.2353514 0.4983295 0.2726721 0.3595736 0.3009704 0.4920244
## d_20 d_21 d_22 d_23 d_25 d_26
## 2008-12-29 0.4446601 0.2684749 0.2532801 0.4174355 0.3996886 0.3506621
## 2009-01-05 0.4439298 0.2795857 0.2700896 0.4237603 0.3972552 0.3569233
## 2009-01-12 0.4764637 0.3266903 0.3268774 0.4580907 0.4298767 0.4033645
## 2009-01-19 0.4271647 0.2844558 0.2849396 0.4231039 0.3764770 0.3515307
## 2009-01-26 0.4140016 0.2809114 0.2861979 0.4183198 0.3612660 0.3435924
## 2009-02-02 0.4083017 0.2847089 0.2960272 0.4194260 0.3541144 0.3451327
## d_28 d_29 d_3 d_31 d_32 d_33
## 2008-12-29 0.2409698 0.1995817 0.3692714 0.5036412 0.1741139 0.2036382
## 2009-01-05 0.2485556 0.2107143 0.3744607 0.5023219 0.1846407 0.2163498
## 2009-01-12 0.2912878 0.2575282 0.4074411 0.5412665 0.2291167 0.2716511
## 2009-01-19 0.2481723 0.2143395 0.3706515 0.4834266 0.1886138 0.2228577
## 2009-01-26 0.2435013 0.2088262 0.3633731 0.4703971 0.1842667 0.2206076
## 2009-02-02 0.2472127 0.2093470 0.3610948 0.4683630 0.1862272 0.2284429
## d_34 d_35 d_36 d_37 d_38 d_39
## 2008-12-29 0.4128579 0.3852355 0.3778185 0.3209316 0.3397267 0.4452858
## 2009-01-05 0.4182664 0.3859322 0.3777464 0.3242407 0.3488059 0.4394305
## 2009-01-12 0.4663383 0.4208596 0.4119223 0.3607388 0.3920069 0.4736202
## 2009-01-19 0.4104876 0.3701521 0.3617756 0.3152357 0.3505981 0.4105697
## 2009-01-26 0.4014505 0.3559547 0.3487988 0.3056072 0.3459112 0.3916376
## 2009-02-02 0.4028516 0.3482007 0.3434432 0.3031027 0.3482623 0.3828116
## d_4 d_40 d_41 d_42 d_43 d_44
## 2008-12-29 0.2337599 0.2196567 0.4394402 0.2474255 0.3687938 0.2565000
## 2009-01-05 0.2593379 0.2403631 0.4318852 0.2683613 0.3694133 0.2686698
## 2009-01-12 0.3592792 0.3351975 0.4576142 0.3636559 0.4042105 0.3208187
## 2009-01-19 0.2814001 0.2517219 0.4015666 0.2811304 0.3545340 0.2742532
## 2009-01-26 0.2866305 0.2506123 0.3817479 0.2817196 0.3417959 0.2709059
## 2009-02-02 0.3128546 0.2695410 0.3695163 0.3033119 0.3363801 0.2761605
## d_45 d_46 d_47 d_48 d_49 d_5
## 2008-12-29 0.4599630 0.2052082 0.2593854 0.4087507 0.3995698 0.3520960
## 2009-01-05 0.4534119 0.2186416 0.2746555 0.4097246 0.4002574 0.3606769
## 2009-01-12 0.4800935 0.2671311 0.3324100 0.4456410 0.4288115 0.4090039
## 2009-01-19 0.4248871 0.2295691 0.2858408 0.3952905 0.3877649 0.3594264
## 2009-01-26 0.4057384 0.2301338 0.2854623 0.3827498 0.3764931 0.3529371
## 2009-02-02 0.3939555 0.2386682 0.2947132 0.3777779 0.3705637 0.3552173
## d_50 d_51 d_52 d_53 d_54 d_55
## 2008-12-29 0.1861633 0.3432858 0.4579974 0.2042064 0.2805190 0.2831990
## 2009-01-05 0.1976868 0.3443546 0.4551332 0.2352984 0.2928957 0.2966881
## 2009-01-12 0.2443557 0.3795868 0.4855001 0.3410345 0.3477630 0.3843874
## 2009-01-19 0.2051665 0.3290406 0.4339684 0.2695353 0.2983259 0.2939665
## 2009-01-26 0.2044018 0.3147844 0.4184742 0.2820350 0.2951119 0.2862090
## 2009-02-02 0.2119916 0.3066881 0.4103175 0.3166709 0.3015707 0.2988774
## d_56 d_57 d_58 d_6 d_60 d_61
## 2008-12-29 0.3236820 0.4069372 0.2873552 0.3110044 0.3325707 0.3335435
## 2009-01-05 0.3325734 0.4028505 0.2997368 0.3172467 0.3407407 0.3441049
## 2009-01-12 0.3812601 0.4347464 0.3480880 0.3636228 0.3914413 0.3906665
## 2009-01-19 0.3321456 0.3773114 0.3070465 0.3116255 0.3379317 0.3479038
## 2009-01-26 0.3262773 0.3585939 0.3046115 0.3033735 0.3307930 0.3438547
## 2009-02-02 0.3293816 0.3472629 0.3094101 0.3044064 0.3335203 0.3471777
## d_62 d_63 d_64 d_65 d_66 d_67
## 2008-12-29 0.2760306 0.3336212 0.4020540 0.3403297 0.2661197 0.2669488
## 2009-01-05 0.2896380 0.3397206 0.3955047 0.3478145 0.2773346 0.2899485
## 2009-01-12 0.3431996 0.3801296 0.4220059 0.3893675 0.3283769 0.3875274
## 2009-01-19 0.2979911 0.3350599 0.3662329 0.3462004 0.2816821 0.3077500
## 2009-01-26 0.2959139 0.3267336 0.3459267 0.3395877 0.2783287 0.3117944
## 2009-02-02 0.3023260 0.3251484 0.3322248 0.3397984 0.2840694 0.3377518
## d_68 d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 0.2688478 0.3396584 0.1683511 NA NA NA NA NA
## 2009-01-05 0.2823005 0.3492668 0.1880529 NA NA NA NA NA
## 2009-01-12 0.3301128 0.3949132 0.2421639 NA NA NA NA NA
## 2009-01-19 0.2925372 0.3513175 0.2109790 NA NA NA NA NA
## 2009-01-26 0.2918368 0.3465619 0.2168041 NA NA NA NA NA
## 2009-02-02 0.2980425 0.3493424 0.2297050 NA NA NA NA NAtail(pred.B.cv$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 1.0150375 1.0776936 1.0198085 1.0591903 0.9841650 0.9620116
## 2020-11-09 NA 0.6251803 0.6913531 0.6313462 0.6727281 0.5863371 0.5937758
## 2020-11-16 NA 1.0243758 1.1019127 1.0535331 1.0886497 1.0044736 1.0171267
## 2020-11-23 NA 0.9527293 1.0288895 0.9916107 1.0192567 0.9329078 0.9577483
## 2020-11-30 NA 1.1483583 1.2294370 1.2083545 1.2296129 1.1446814 1.1793494
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2020-11-02 1.1561735 1.147123 1.164793 1.020855 1.116269 1.2825994 1.0213387
## 2020-11-09 0.7553162 0.723739 0.758770 0.643016 0.725545 0.8196408 0.6443358
## 2020-11-16 1.1506762 1.125905 1.152468 1.067550 1.125433 1.1909961 1.0745858
## 2020-11-23 1.0582438 1.017714 1.065531 1.007662 1.043340 1.0731742 1.0132734
## 2020-11-30 1.2396615 1.214072 1.253733 1.224595 1.228521 1.2290675 1.2376935
## 2020-12-07 NA NA NA NA NA NA NA
## d_21 d_22 d_23 d_25 d_26 d_28
## 2020-11-02 1.0448420 1.0920563 1.0911951 0.9647147 1.0342517 0.9003043
## 2020-11-09 0.6539470 0.6925186 0.6945855 0.5946751 0.6377978 0.5131376
## 2020-11-16 1.0615172 1.0932745 1.1080573 1.0242741 1.0620301 0.9198463
## 2020-11-23 0.9891587 1.0159927 1.0238217 0.9666114 1.0002494 0.8438110
## 2020-11-30 1.1895191 1.2100321 1.2224294 1.1876949 1.2141364 1.0445219
## 2020-12-07 NA NA NA NA NA NA
## d_29 d_3 d_31 d_32 d_33 d_34
## 2020-11-02 1.0361800 1.0783991 1.0059734 0.9684093 0.9456605 1.0617426
## 2020-11-09 0.6552218 0.6833334 0.6340943 0.5773669 0.5591726 0.6839008
## 2020-11-16 1.0648836 1.0983172 1.0610828 0.9879177 0.9621700 1.1041349
## 2020-11-23 0.9970938 1.0236954 1.0015543 0.9135634 0.8921258 1.0434900
## 2020-11-30 1.2023364 1.2309221 1.2138238 1.1171841 1.0893082 1.2555070
## 2020-12-07 NA NA NA NA NA NA
## d_35 d_36 d_37 d_38 d_39 d_4
## 2020-11-02 1.0620075 0.9817250 0.9780876 1.0803897 0.9236184 1.0462908
## 2020-11-09 0.6892449 0.6046740 0.5928632 0.6885219 0.5715089 0.6454975
## 2020-11-16 1.1178506 1.0294449 1.0136770 1.0997828 0.9983757 1.0551971
## 2020-11-23 1.0615771 0.9690825 0.9455575 1.0237380 0.9552188 0.9971032
## 2020-11-30 1.2839672 1.1857704 1.1568365 1.2283501 1.1776324 1.1925871
## 2020-12-07 NA NA NA NA NA NA
## d_40 d_41 d_42 d_43 d_44 d_45
## 2020-11-02 0.9054965 0.9479669 1.0124860 0.9970620 1.0414596 0.9870247
## 2020-11-09 0.5470347 0.5803466 0.6255479 0.6206986 0.6487289 0.6205416
## 2020-11-16 0.9416465 1.0204236 1.0446210 1.0461720 1.0589809 1.0570829
## 2020-11-23 0.8801214 0.9676142 0.9917533 0.9838337 0.9898298 1.0052004
## 2020-11-30 1.0765392 1.1988128 1.1928202 1.2017317 1.1926276 1.2333951
## 2020-12-07 NA NA NA NA NA NA
## d_46 d_47 d_48 d_49 d_5 d_50
## 2020-11-02 0.9517694 1.0548802 1.0339263 1.0387043 1.0456600 0.8876971
## 2020-11-09 0.5529526 0.6716038 0.6595257 0.6528972 0.6698027 0.4947756
## 2020-11-16 0.9489201 1.0755656 1.0816021 1.0767264 1.0721795 0.8906770
## 2020-11-23 0.8679939 1.0014233 1.0234128 1.0070041 1.0084142 0.8111266
## 2020-11-30 1.0595220 1.1872085 1.2388418 1.2182878 1.2118826 1.0028663
## 2020-12-07 NA NA NA NA NA NA
## d_51 d_52 d_53 d_54 d_55 d_56
## 2020-11-02 1.0327922 1.0314674 1.1132033 1.0377687 0.9845357 1.0231886
## 2020-11-09 0.6669587 0.6537372 0.6812753 0.6531501 0.6206208 0.6420580
## 2020-11-16 1.0959861 1.0842459 1.0826501 1.0651263 1.0527421 1.0461825
## 2020-11-23 1.0403651 1.0286754 1.0082953 0.9966882 1.0133027 0.9831596
## 2020-11-30 1.2553539 1.2509919 1.1861315 1.2032899 1.2285154 1.1868197
## 2020-12-07 NA NA NA NA NA NA
## d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 0.9935330 1.0948235 0.9957890 1.0276878 1.0841772 1.0625656
## 2020-11-09 0.6373275 0.6998669 0.6120707 0.6484534 0.6950927 0.6726742
## 2020-11-16 1.0673926 1.1071086 1.0330153 1.0615123 1.0978030 1.0729691
## 2020-11-23 1.0219978 1.0326973 0.9690720 0.9982060 1.0288889 1.0037612
## 2020-11-30 1.2476433 1.2302258 1.1791473 1.2061263 1.2307727 1.2029491
## 2020-12-07 NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68
## 2020-11-02 1.0552428 0.9614875 1.0744089 1.0073493 1.0145668 1.0898587
## 2020-11-09 0.6727983 0.6009151 0.6723145 0.6174771 0.6280493 0.6869631
## 2020-11-16 1.0878553 1.0364567 1.0906425 1.0263662 1.0443520 1.0910055
## 2020-11-23 1.0226460 0.9897999 1.0246298 0.9530366 0.9796384 1.0119553
## 2020-11-30 1.2316994 1.2194697 1.2377665 1.1550219 1.1658121 1.2070013
## 2020-12-07 NA NA NA NA NA NA
## d_7 d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 1.0722112 1.0336011 NA NA NA NA NA
## 2020-11-09 0.6829798 0.6237523 NA NA NA NA NA
## 2020-11-16 1.0884182 1.0135626 NA NA NA NA NA
## 2020-11-23 1.0188960 0.9224215 NA NA NA NA NA
## 2020-11-30 1.2223381 1.1060283 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NA NAhead(pred.B.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2008-12-29 NA 1.282235 1.442799 1.583269 1.539445 1.395196 1.479022
## 2009-01-05 NA 1.300224 1.456188 1.583662 1.542359 1.399704 1.473440
## 2009-01-12 NA 1.371921 1.523724 1.637776 1.588965 1.443938 1.509156
## 2009-01-19 NA 1.312506 1.458428 1.560892 1.526900 1.389425 1.441364
## 2009-01-26 NA 1.310834 1.451441 1.542684 1.511546 1.377356 1.417528
## 2009-02-02 NA 1.320527 1.455698 1.536513 1.504268 1.372661 1.400658
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2008-12-29 1.213964 1.626691 1.271460 1.519599 1.313545 1.476867 1.602453
## 2009-01-05 1.241845 1.646174 1.297752 1.520258 1.339146 1.522727 1.600861
## 2009-01-12 1.314912 1.713899 1.370525 1.573836 1.424383 1.635405 1.653416
## 2009-01-19 1.275833 1.666266 1.328104 1.497353 1.364677 1.594005 1.573571
## 2009-01-26 1.283467 1.671668 1.334058 1.478216 1.368828 1.626395 1.552746
## 2009-02-02 1.298132 1.689626 1.347634 1.470061 1.385670 1.679228 1.543732
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2008-12-29 1.342389 1.323985 1.559910 1.530633 1.459404 1.305992 1.253021
## 2009-01-05 1.357173 1.346040 1.569373 1.526742 1.468146 1.315789 1.266849
## 2009-01-12 1.422434 1.424330 1.623792 1.577212 1.537547 1.373098 1.327381
## 2009-01-19 1.363448 1.365540 1.567641 1.495072 1.459569 1.315044 1.271122
## 2009-01-26 1.358491 1.367026 1.559902 1.472401 1.447781 1.308825 1.264010
## 2009-02-02 1.363556 1.380350 1.561433 1.461831 1.449824 1.313622 1.264573
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2008-12-29 1.486559 1.699155 1.222324 1.258190 1.550991 1.509530 1.498491
## 2009-01-05 1.493880 1.696522 1.234951 1.274004 1.559057 1.510217 1.498009
## 2009-01-12 1.543595 1.763537 1.290835 1.346180 1.635516 1.563566 1.549752
## 2009-01-19 1.487534 1.664141 1.239369 1.281859 1.546419 1.485993 1.473682
## 2009-01-26 1.476503 1.642368 1.233811 1.278805 1.532299 1.464834 1.454468
## 2009-02-02 1.472958 1.638853 1.236095 1.288729 1.534286 1.453361 1.446538
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2008-12-29 1.415959 1.442484 1.602837 1.297257 1.279090 1.594110 1.315106
## 2009-01-05 1.420277 1.455278 1.593111 1.330500 1.305492 1.581694 1.342560
## 2009-01-12 1.472731 1.519193 1.648186 1.469989 1.435012 1.622541 1.476437
## 2009-01-19 1.406942 1.457302 1.547208 1.359575 1.319818 1.533802 1.359208
## 2009-01-26 1.393238 1.450274 1.517979 1.366480 1.318138 1.503463 1.359786
## 2009-02-02 1.389583 1.453526 1.504477 1.402613 1.343157 1.485004 1.389292
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2008-12-29 1.485028 1.328255 1.627164 1.260113 1.330323 1.544567 1.532289
## 2009-01-05 1.485578 1.344130 1.616112 1.277011 1.350494 1.545898 1.532918
## 2009-01-12 1.537846 1.415726 1.659428 1.340326 1.430508 1.602272 1.576937
## 2009-01-19 1.463049 1.351026 1.569991 1.290806 1.365192 1.523466 1.513210
## 2009-01-26 1.444318 1.346280 1.539968 1.291446 1.364490 1.504376 1.496002
## 2009-02-02 1.436357 1.353197 1.521743 1.302445 1.377027 1.496835 1.486972
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2008-12-29 1.460218 1.236341 1.447515 1.623969 1.259479 1.358736 1.363004
## 2009-01-05 1.472461 1.250530 1.448713 1.618896 1.298897 1.375353 1.381134
## 2009-01-12 1.545054 1.310145 1.500347 1.668424 1.443413 1.452640 1.507365
## 2009-01-19 1.470068 1.259688 1.426140 1.584314 1.343539 1.382343 1.376768
## 2009-01-26 1.460354 1.258638 1.405751 1.559706 1.360215 1.377721 1.365904
## 2009-02-02 1.463530 1.268160 1.394264 1.546847 1.407976 1.386502 1.383145
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2008-12-29 1.419311 1.541743 1.367975 1.402660 1.431334 1.432641 1.354452
## 2009-01-05 1.431656 1.535212 1.384799 1.411036 1.442756 1.447623 1.372612
## 2009-01-12 1.502778 1.584749 1.453199 1.477641 1.517495 1.516410 1.447769
## 2009-01-19 1.430506 1.496115 1.394599 1.402472 1.438191 1.452759 1.383480
## 2009-01-26 1.421936 1.468229 1.391072 1.390708 1.427767 1.446747 1.380372
## 2009-02-02 1.426203 1.451576 1.397657 1.391965 1.431516 1.451453 1.389072
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2008-12-29 1.433593 1.535613 1.443354 1.339919 1.341034 1.343678 1.441428
## 2009-01-05 1.442016 1.525186 1.453836 1.354718 1.371857 1.361547 1.455115
## 2009-01-12 1.501162 1.565782 1.515189 1.425371 1.512106 1.427925 1.522864
## 2009-01-19 1.434750 1.480553 1.450907 1.360109 1.395876 1.375020 1.457727
## 2009-01-26 1.422649 1.450561 1.441132 1.355365 1.401302 1.373860 1.450669
## 2009-02-02 1.420241 1.430646 1.441276 1.363020 1.437972 1.382261 1.454598
## d_8 d_24 d_27 d_30 d_59 d_9
## 2008-12-29 1.215206 NA NA NA NA NA
## 2009-01-05 1.239086 NA NA NA NA NA
## 2009-01-12 1.307705 NA NA NA NA NA
## 2009-01-19 1.267327 NA NA NA NA NA
## 2009-01-26 1.274548 NA NA NA NA NA
## 2009-02-02 1.290957 NA NA NA NA NAtail(pred.B.cv.log$pred.all$EX)## central d_1 d_10 d_11 d_12 d_13 d_14
## 2020-11-02 NA 2.774707 2.950054 2.785553 2.898576 2.687846 2.630276
## 2020-11-09 NA 1.878626 2.004501 1.888545 1.969049 1.805260 1.819653
## 2020-11-16 NA 2.800152 3.022040 2.880341 2.984271 2.742101 2.778414
## 2020-11-23 NA 2.606599 2.809342 2.707406 2.784120 2.552640 2.618159
## 2020-11-30 NA 3.170127 3.433569 3.362966 3.436103 3.154885 3.267818
## 2020-12-07 NA NA NA NA NA NA NA
## d_15 d_16 d_17 d_18 d_19 d_2 d_20
## 2020-11-02 3.190019 3.160998 3.220468 2.788032 3.065765 3.621659 2.788091
## 2020-11-09 2.136362 2.069476 2.145388 1.910367 2.073920 2.279030 1.912050
## 2020-11-16 3.172268 3.093627 3.180104 2.920361 3.093432 3.303470 2.939799
## 2020-11-23 2.892275 2.776311 2.915228 2.750464 2.849728 2.936171 2.764980
## 2020-11-30 3.467849 3.378837 3.519052 3.416862 3.429857 3.431631 3.460948
## 2020-12-07 NA NA NA NA NA NA NA
## d_21 d_22 d_23 d_25 d_26 d_28 d_29
## 2020-11-02 2.855439 2.997261 2.989548 2.635053 2.826848 2.471497 2.831633
## 2020-11-09 1.931396 2.009606 2.010341 1.819930 1.901200 1.677983 1.934420
## 2020-11-16 2.903134 2.999853 3.039397 2.796531 2.905443 2.520062 2.913769
## 2020-11-23 2.700599 2.776625 2.793773 2.639914 2.731252 2.335620 2.722903
## 2020-11-30 3.300045 3.371342 3.407747 3.293342 3.382729 2.854969 3.343584
## 2020-12-07 NA NA NA NA NA NA NA
## d_3 d_31 d_32 d_33 d_34 d_35 d_36
## 2020-11-02 2.953556 2.746222 2.645493 2.587066 2.903588 2.905010 2.680959
## 2020-11-09 1.989208 1.893068 1.788923 1.757522 1.989679 2.000687 1.838449
## 2020-11-16 3.011998 2.901200 2.696736 2.629663 3.028781 3.070965 2.811085
## 2020-11-23 2.795333 2.733588 2.503372 2.451859 2.850661 2.902839 2.646291
## 2020-11-30 3.439163 3.380359 3.068778 2.986619 3.524298 3.625967 3.286644
## 2020-12-07 NA NA NA NA NA NA NA
## d_37 d_38 d_39 d_4 d_40 d_41 d_42
## 2020-11-02 2.670951 2.958414 2.532819 2.862185 2.485235 2.591327 2.766262
## 2020-11-09 1.816719 1.998874 1.780825 1.916678 1.736214 1.793864 1.878288
## 2020-11-16 2.766993 3.015369 2.728852 2.886909 2.575932 2.785324 2.855727
## 2020-11-23 2.584795 2.794431 2.613638 2.723922 2.422177 2.642061 2.708619
## 2020-11-30 3.193170 3.428972 3.264977 3.312219 2.948057 3.329583 3.312052
## 2020-12-07 NA NA NA NA NA NA NA
## d_43 d_44 d_45 d_46 d_47 d_48 d_49
## 2020-11-02 2.722106 2.849035 2.695286 2.602890 2.883439 2.825016 2.837492
## 2020-11-09 1.867950 1.923271 1.867953 1.746722 1.965163 1.942641 1.928815
## 2020-11-16 2.858200 2.898368 2.890126 2.595271 2.943178 2.962719 2.946500
## 2020-11-23 2.685333 2.704593 2.744014 2.393588 2.732954 2.795321 2.747977
## 2020-11-30 3.339171 3.312770 3.447694 2.899082 3.291298 3.467551 3.394645
## 2020-12-07 NA NA NA NA NA NA NA
## d_5 d_50 d_51 d_52 d_53 d_54 d_55
## 2020-11-02 2.861379 2.441746 2.820986 2.817935 3.056913 2.834571 2.690471
## 2020-11-09 1.964621 1.648271 1.956328 1.931113 1.984325 1.929266 1.869373
## 2020-11-16 2.937661 2.448831 3.004143 2.969877 2.964019 2.912666 2.879503
## 2020-11-23 2.756239 2.261637 2.841529 2.809354 2.751570 2.720093 2.768092
## 2020-11-30 3.378500 2.739844 3.523211 3.509101 3.287312 3.344715 3.432994
## 2020-12-07 NA NA NA NA NA NA NA
## d_56 d_57 d_58 d_6 d_60 d_61 d_62
## 2020-11-02 2.798000 2.714004 3.000952 2.720691 2.806878 2.971894 2.911998
## 2020-11-09 1.911002 1.900529 2.021591 1.853257 1.920733 2.013808 1.971365
## 2020-11-16 2.862483 2.921727 3.037711 2.822880 2.902924 3.012335 2.941412
## 2020-11-23 2.687699 2.792176 2.819993 2.647907 2.724943 2.811850 2.744603
## 2020-11-30 3.295118 3.499320 3.436196 3.267027 3.355099 3.441227 3.349662
## 2020-12-07 NA NA NA NA NA NA NA
## d_63 d_64 d_65 d_66 d_67 d_68 d_7
## 2020-11-02 2.886553 2.628531 2.941277 2.750721 2.769566 2.986962 2.935527
## 2020-11-09 1.968822 1.832489 1.967076 1.862357 1.881321 1.996064 1.988873
## 2020-11-16 2.981386 2.832422 2.988446 2.802941 2.852425 2.989528 2.983163
## 2020-11-23 2.793093 2.703318 2.797404 2.604808 2.673627 2.762224 2.782926
## 2020-11-30 3.442660 3.401574 3.462007 3.188151 3.220937 3.357255 3.411142
## 2020-12-07 NA NA NA NA NA NA NA
## d_8 d_24 d_27 d_30 d_59 d_9
## 2020-11-02 2.824528 NA NA NA NA NA
## 2020-11-09 1.874438 NA NA NA NA NA
## 2020-11-16 2.767696 NA NA NA NA NA
## 2020-11-23 2.526525 NA NA NA NA NA
## 2020-11-30 3.035855 NA NA NA NA NA
## 2020-12-07 NA NA NA NA NA NAsummary(pred.B.cv)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX
## obs 0.13784779 0.13784779 0.11072223
## average 0.08626559 0.08626559 0.06637898
##
## R2:
## EX.mu EX.mu.beta EX
## obs 0.7179091 0.7179091 0.8180051
## average 0.6694103 0.6694103 0.8042619
##
## Coverage of 95% prediction intervals:
## EX
## obs 0.9367089
## average 0.8571429summary(pred.B.cv.log)## Cross-validation predictions for STmodel with 10 CV-groups.
## Predictions for 790 observations.
## Temporal averages for 63 locations.
##
## RMSE:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.1938374 0.1938374 0.14992738 0.14995756
## average 0.1017864 0.1017864 0.06827081 0.06883603
##
## R2:
## EX.mu EX.mu.beta EX EX.pred
## obs 0.6931546 0.6931546 0.8164280 0.8163541
## average 0.6921204 0.6921204 0.8614931 0.8591901
##
## Coverage of 95% prediction intervals:
## EX.pred
## obs 0.9367089
## average 0.8730159par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
jpeg(filename = here::here("Figs", "ST_CV_Obs_vs_Pred_BC_ModB.jpeg"),
width = 8, height = 4, units = "in", res = 500)
par(mfrow=c(1,2), mar=c(3.3,3.3,1.5,1), mgp=c(2,1,0))
plot(pred.B.cv, "obs", ID="all", pch=c(19,NA), cex=.25, lty=c(NA,2),
col=c("ID", "black", "grey"),
ylim=c(-1,2),
xlab="Observations", ylab="Predictions",
main="Cross-validation BC (log ug/m3)")
with(pred.B.cv.log$pred.LTA, plotCI(obs, EX.pred, uiw=1.96*sqrt(VX.pred),
xlab="Observations", ylab="Predictions",
main="Temporal average BC (ug/m3)"))
abline(0, 1, col="grey")
dev.off()## quartz_off_screen
## 23.2.5 Predictions for 2009-2020
What do the long-term predictions look like for this model? Predicting at the distributed (residential + community) sites
x.B <- coef(est.denver.model.B, pars = "cov")$par
x.B## [1] -10.463290 -7.028546 -7.730966 13.512563 -2.978499 -4.934033E.B <- predict(denver.model.B, est.denver.model.B, LTA = T,
transform="unbiased", pred.var=FALSE)
print(E.B)## Prediction for STmodel.
##
## Regression parameters:
## 0 Spatio-temporal covariate(s).
## 42 beta-fields regression parameters in x$pars.
##
## Regression parameters are assumed to be known and
## prediction variances do NOT include
## uncertainties in regression parameters.
##
## Prediction of beta-fields, (x$beta):
## List of 2
## $ mu: num [1:69, 1:3] 0.403 0.203 0.276 0.237 0.252 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX: num [1:69, 1:3] 0.404 0.204 0.275 0.236 0.25 ...
## ..- attr(*, "dimnames")=List of 2
##
## Predictions for log-Gaussian field of type: unbiased
##
## Predictions for 622 times at 69 locations.
## List of 5
## $ EX.mu : num [1:622, 1:69] 1.5 1.53 1.62 1.56 1.58 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.mu.beta: num [1:622, 1:69] 1.54 1.57 1.65 1.6 1.61 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX : num [1:622, 1:69] 1.58 1.61 1.7 1.64 1.65 ...
## ..- attr(*, "dimnames")=List of 2
## $ EX.pred : num [1:622, 1:69] 1.59 1.61 1.7 1.64 1.65 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.EX : num [1:622, 1:69] 0.432 0.448 0.503 0.467 0.474 ...
## ..- attr(*, "dimnames")=List of 2
##
## Variances have been computed.
## List of 2
## $ log.VX : num [1:622, 1:69] 0.0514 0.0513 0.0512 0.0511 0.0511 ...
## ..- attr(*, "dimnames")=List of 2
## $ log.VX.pred: num [1:622, 1:69] 0.0586 0.0585 0.0584 0.0583 0.0583 ...
## ..- attr(*, "dimnames")=List of 2
##
## Mean squared prediciton errors NOT been computed.
##
## 69 temporal averages have been compute.
## List of 1
## $ LTA:'data.frame': 69 obs. of 4 variables:week_preds <- as.data.frame(E.B$EX) %>%
mutate(week = as.Date(rownames(E.B$EX)),
year = year(as.Date(rownames(E.B$EX))))
week_sites <- colnames(week_preds)[which(str_detect(colnames(week_preds), "d_"))]
box_preds <- week_preds %>%
select(week, all_of(week_sites)) %>%
#filter(week %in% week_weeks) %>%
mutate(month = month(week),
year = year(week)) %>%
filter(year %in% c(2009:2020))
box_data <- box_preds %>%
pivot_longer(-c(week, month, year), names_to = "location", values_to = "pred")
summary(box_data)## week month year location
## Min. :2009-01-05 Min. : 1.000 Min. :2009 Length:42228
## 1st Qu.:2012-01-09 1st Qu.: 4.000 1st Qu.:2012 Class :character
## Median :2014-12-29 Median : 7.000 Median :2014 Mode :character
## Mean :2014-12-28 Mean : 6.494 Mean :2014
## 3rd Qu.:2017-12-18 3rd Qu.: 9.000 3rd Qu.:2017
## Max. :2020-12-07 Max. :12.000 Max. :2020
##
## pred
## Min. :0.3553
## 1st Qu.:1.1765
## Median :1.4744
## Mean :1.5943
## 3rd Qu.:2.0065
## Max. :4.0346
## NA's :861#' Box plot of weekly BC predicted at all distributed sites grouped by month
box_summary <- box_data %>%
group_by(month) %>%
summarize(median = median(pred, na.rm=T)) %>%
arrange(desc(median))
box_summarypred_box_plot <- ggplot(box_data) +
geom_boxplot(aes(x = as.factor(month), y = pred), #, color = as.factor(month)),
show.legend = F) +
#scale_color_viridis(name = "Month", discrete = T) +
facet_wrap(. ~ year, scales = "free") +
xlab("") + ylab("Distributed site BC concentration (\u03BCg/m\u00B3): Predicted") +
# theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
scale_x_discrete(labels = str_sub(month.name, 1, 1))
pred_box_plot
4 Summary Table
| model | trend_data | cov_beta0 | cov_beta1 | cov_beta2 | cov_nu_error | no_basis_fns | df_per_year | all_converge | cv_rmse_obs | cv_rmse_avg | cv_r2_obs | cv_r2_avg | cv_coverage_obs | cv_coverage_avg |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Model A | BC + NO2 + PM2.5 | iid | iid | NA | exp | 1 | 4 | TRUE | 0.15 | 0.06 | 0.81 | 0.89 | 0.94 | 0.90 |
| Model B | BC + NO2 + PM2.5 | iid | iid | iid | exp | 2 | 4 | TRUE | 0.15 | 0.07 | 0.82 | 0.86 | 0.94 | 0.87 |